Pharmacophore Modeling for Pleurotin Derivatives Targeting Human Thioredoxin Reductase

Abstract

Pleurotin is a natural compound and potent inhibitor of thioredoxin reductase (TrxR) enzyme, which is a target for cancer chemotherapy. The present study proposes a pharmacophore model for TrxR inhibitors based on the in silico analysis of 34 pleurotin analogues. The results led to ligand 33 with a flexible -SR side-chain at p-quinone and a H-bond donor -CH₂NH₂ at the oxepane ring. The ligand-TrxR binding free energy was -236 kJ mol 1 for 33, better than pleurotin (-166 kJ mol 1). The results provide a clear and practical guide to design pleurotin derivatives with anticancer potential that could expand the pharmacological potential of this class of natural products.

Keywords:
pleurotin; TrxR; pharmacophore model; molecular docking; molecular dynamics

Introduction

According to the International Agency for Research on Cancer (IARC),¹ in 2020, ca.19 million people were diagnosed with cancer in the world, with ca. 52% deaths.¹ Currently, there is a range of possibilities for cancer treatment, among them, are targeted therapy (chemotherapy), radiation therapy, immunotherapy, and surgery. Despite the advances in medicine, the mortality remains high and the side effects are marked for most treatments. In view of this, targeted therapy is a promising way to treat cancer with reduced side effects and high specificity.² This strategy mostly consists of using small molecules to inhibit biomolecules (proteins and nucleotides) that are associated with cell proliferation and survival, or use monoclonal antibodies to help the immune system to fight against cancer.

At this point, the thioredoxin system, which comprehends thioredoxin (Trx), thioredoxin reductase (TrxR) and the cofactor nicotinamide adenine dinucleotide phosphate (NADPH) was associated with high expressiveness in tumor cells and is a potential target for small inhibitors.³,⁴,⁵,⁶ The TrxR is part of the cell antioxidant system and is found in humans in three main forms, cytoplasmic (TrxR1), mitochondrial (TrxR2) and a testis-specific thioredoxin glutathione reductase (TGR or TrxR3).⁷ The present study focuses on TrxR1 (hereinafter referred to as TrxR), which is characterized to catalyze the reduction of Trx, and other oxidized substrates, through the transport of electrons from NADPH to the substrate, promoting the breaking and forming of disulfide bonds between Cys in the N- and C-terminal moieties of TrxR (Figure 1).⁷ The TrxR binding site is described as a cavity ca. 15 Å deep on the enzyme surface with an open “mouth” to the solvent, from which the inhibitors can access and interact with the enzyme.

Figure 1
Structure of TrxR (PDB code: 2J3N) with emphasis on the binding site. The enzyme is composed by two chains, highlighted in blue (chain D) and green (chain C), and two cofactors, FADH₂ (red) and NADP⁺ (yellow) each. In the N-terminal, Cys59 and Cys64 are in an oxidized state, forming a disulfide bond. In the C-terminal, Cys497 and Cys498 are in the reduced state, in thiol form.

Recent studies about TrxR inhibition suggest two broad groups of potential molecules; the metal complexes and the organic compounds. The metal complexes include molecules based on transition metals, mainly Au, like auranofin (half maximal inhibitory concentration (IC₅₀) ca.1.5 μM), and Gd, like motexafin gadolinium (IC₅₀: 6 μM for rat TrxR1), and others, that were expected to interact covalently with the C-terminal moiety.⁸,⁹,¹⁰ Within the organic inhibitors, natural-derived compounds like curcumins (IC₅₀: 0.3-62 μM),¹¹ benzotriazinones (0.23-5 μM) and naphthoquinone spiroketals (0.27-3.2 μM) are promising TrxR inhibitors, which are expected to interact with the TrxR by forming a non-covalent network of contacts with the catalytic site, blocking the communication of N- and C-terminal domains.¹²,¹³ Focusing on natural-derived compounds, despite the poor structural variability, the molecule named pleurotin stands out as a potential inhibitor of the TrxR-Trx system (IC₅₀: 0.17 μM).¹⁴,¹⁵ This compound was identified in our previous paper¹⁶ as the best representative among 72 molecules evaluated in silico as TrxR inhibitors. In addition to support the potential of pleurotin as TrxR inhibitor, we characterized the binding mode, which involved a large network of drug-receptor contacts with the residues in the cavity, in particular, Tyr116 and His472, both being essential for enzyme catalytic function.¹⁷ Recently, Er-rajy et al.,¹⁸ proposed four new natural-derived inhibitors for TrxR based on quantitative structure-activity relationship (QSAR), docking, and molecular dynamics (MD) analysis; however, different from pleurotin binding mode,¹⁶ the authors proposed the flavin adenine dinucleotide (FAD) domain as binding site.

In the present study, we extend the investigation of pleurotin as inhibitor of TrxR aiming to construct a pleurotin-based pharmacophore. A set of 14 known natural analogues of pleurotin were evaluated through docking and MD simulations, then, 20 new analogues were designed and evaluated in silico as TrxR inhibitors.

Methodology

Ligand’s collection

A set of 14 pleurotin derivatives were selected from Sandargo et al.¹⁹,²⁰ (Figure 2). The 3D structures were constructed and optimized as a minimum point on the potential energy surface (PES) in gas phase at the semi-empirical level Austin Model 1 (AM1), using the Gaussian 03 program.²¹,²² In addition to the initial set of compounds, 20 new analogues were proposed based on the established pharmacophore, and their geometries optimized at the very same level of theory applied for pleurorin (Figure 3).

Figure 2
Structure of pleurotin (1) and its known natural analogues: leucopleurotin (2), dihydropleurotinic acid (3), pleurogrisein (4), 3-hydroxy-dihydropleurotinic acid (5), 14-hydroxy-dihydropleurotinic acid (6), leucopleurotinic acid (7), 14-oxo-leucopleurotinic acid (8), nematoctone (9), di-oxo-leucopleurotinic acid (10), 4-hydroxy-pleurogrisein (11), thiopleurotinic acid A (12), thiopleurotinic acid B (13) and pleurothiazole (14).

Figure 3
Structure of the 20 new pleurotin analogues proposed in the present work based on the pharmacophore model.

Receptor preparation

The crystal structure of human TrxR was downloaded from Protein Data Bank²³ (PDB, code: 2J3N) and edited using the Visual Molecular Dynamics (VMD)²⁴ and Pymol programs.²⁵ Only the chains C and D and their cofactors were selected (see Figure 1). The FADH₂ was in the reduced form, and the NADP⁺ was in the oxidized state. The amino acids were in their standard forms at neutral pH, with exception of Cys at the N-terminal that were in the oxidized state forming a disulfide bond. After that, the structure of TrxR was deprotonated and renumbered in ascending order (chain C: Ser1-Gly490; chain D: Lys493-Gly983), with pdb4amber utility of the AMBER 16 package.²⁶ The original numbering sequence (PDB code: 2J3N) is Ser10-Gly499 for chain C and Lys9-Gly499 for chain D.

Molecular docking

Docking simulations were carried out for all compounds shown in Figures 2 and 3 using the Autodock 4.2.6.²⁷ The receptor was protonated and the docking space was defined as a 40 Å edge cube centered at the catalytic site of TrxR, i.e., at the S of Cys489 at the C-terminal arm (x, y, z: -22.444, 30.700, 30.894). The docking parameters were kept in their default values, and the best 10 poses were selected for each compound according to the ligand-receptor binding energy (B.E.).

Molecular dynamics (MD) simulation

To evaluate the stability and the ligand-receptor interactions under thermodynamic conditions (aqueous solution, T = 310 K and p = 1 atm), MD simulations were performed for the best pose in docking. These calculations were performed using the AMBER 16 program with the ff14SB force-field for the receptor and the GAFF2 (AM1-BCC atomic charges) for the ligands and cofactors.²⁸,²⁹,³⁰ The solvent was represented for the Self-Consistent Reaction Field (SCRF) continuum model with dielectric constant of 78.3553 (water).

The simulation protocol was: (i) optimization of the entire system, (ii) staggered heating from 50 to 310 K for 3 ns (τ = 2 fs), with the ligand restricted with a force constant of 10.0 kcal mol^-1 Å^-2, (iii) equilibration at 310 K for 1.4 ns (τ = 1 fs), with decreasing ligand constraint to zero, (iv) two production runs of 25 ns each (τ = 2 fs). The temperature was controlled using a Langevin thermostat.³¹ The H-bonds were constrained with the SHAKE algorithm and the non-bonded cutoff was 9999.9Å (effectively infinite). The last 25 ns of the production trajectories were analyzed using the cpptraj utility, and were employed to predict the ligand-receptor binding free energy (AG_aq) using the Generalized Born/Surface Area model (GB/SA) with option igb = 2.

Results and Discussion

Preliminary analyses

The pleurotin (1) was taken as a template to stablish the minimum pharmacophore based on the analysis of the binding mode with TrxR.¹⁶ We observed the formation of a broad and strong network of non-covalent contacts (Figure 4), which keep the ligand stable between the C- and the N-terminal motifs, blocking the movement of the C-terminal domain towards the N-terminal that is essential for enzyme catalytic function. Within this network of contacts, the hot-spots residues His463, Tyr600, Cys489, Ile831 and Glu468 are in short contact with the ligand and, therefore, are responsible to anchor the ligand at the binding site.

Figure 4
Network of residues which form non-covalent contacts with the pleurotin (1).¹⁶ The highlighted residues were assigned as hot-spot residues (residues with a B.E. contribution lower than -4 kJ mol^-1) or have a frequency of contacts with the ligand > 50% in the MD trajectory. The 2D (left) and 3D (right) structures are shown.

In order to clarify the role of different molecular groups on the pleurotin structure for the binding mode with TrxR, a series of in silico tests were carried out for 14 known pleurotin analogues shown in Figure 2. The docking simulations were first applied. The results showed three compounds with high score (B.E. < –31 kJ mol^-1 in the last column of Table 1), 2 (leucopleurotin; B.E = –35.57 kJ mol^-1), 9 (nematoctone; B.E = –35.57 kJ mol^-1) and the pleurotin itself 1 (B.E = –40.09 kJ mol^-1) (see also Figure S1 in the Supplementary Information (SI) section). The other structures presented B.E. within the range of –29 ± 2 kJ mol^-1. Figure 2 shows that the selected pleurotin derivatives are very similar, except the compounds 12, 13 and 14, which have an -SR group in their structures at G4 position. In view of the structural similarities, docking results alone may not be sufficient to provide a good understanding of the necessary features to the binding mode and screening of derivatives. Thus, MD simulations were conducted for all 14 derivatives using the best poses from docking as initial guess structures.

Table 1
Parameters calculated from docking and MD trajectories for the 14 compounds shown in Figure 2

	MD simulation					Docking
	RMSD / Å	dCN / Å	dCL / Å	dNL / Å	ΔGaq (GB/SA) / (kJ mol1)	B.E. / (kJ mol-1)
Pleurotin 1	1.6 ± 0.4	14 ± 1	5 ± 1	9.7 ± 0.9	-166 ± 21 (±1.3)	-40.09
Analogues of leucopleurotin
2	1.4 ± 0.2	11 ± 1	9.5 ± 0.9	12.4 ± 0.5	-169 ± 14 (±0.9)	-35.57
7	1.7 ± 0.3	17 ± 1	8 ± 2	9.9 ± 0.7	-133 ± 24 (±1.5)	-28.05
8	1.4 ± 0.2	19 ± 2	9 ± 3	10.9 ± 0.7	-149 ± 16 (±0.2)	-27.96
9	1.4 ± 0.2	17 ± 1	10 ± 1	8.9 ± 0.5	-143 ± 14 (±0.9)	-35.57
Analogues of dihydropleurotinic acid
3	1.7 ± 0.3	22.4 ± 0.9	11.1 ± 0.5	12.1 ± 0.8	-159 ± 15 (±1.0)	-30.81
5	1.5 ± 0.2	16.5 ± 0.9	8 ± 1	8.3 ± 0.7	-158 ± 19 (±1.0)	-30.14
6	1.7 ± 0.2	12 ± 1	6.7 ± 0.6	12 ± 1	-154 ± 16 (±1.0)	-27.88
Analogues of pleurogrisein
4	1.5 ± 0.2	14 ± 3	8 ± 2	10.6 ± 0.7	-149 ± 20 (±1.3)	-29.76
10	1.3 ± 0.2	17.2 ± 0.9	8 ± 1	11.8 ± 0.5	-163 ± 18 (±1.1)	-25.16
11	1.4 ± 0.2	12.9 ± 0.9	7.2 ± 0.8	9.3 ± 0.5	-128 ± 13 (±0.8)	-28.51
S-bearing analogues
12	1.5 ± 0.2	15.9 ± 0.9	8.8 ± 0.9	14 ± 2	-155 ± 21 (±1.3)	-29.64
13	1.3 ± 0.1	14.1 ± 0.6	3.9 ± 0.4	11.2 ± 0.5	-211 ± 17 (±1.1)	-29.26
14	1.4 ± 0.2	21 ± 1	10.5 ± 0.8	11.3 ± 0.7	-155 ± 13 (±0.8)	-29.05

• The root mean square deviation (RMSD) refers to the protein fluctuation only. The ΔG_aq (= ΔG_v.d.w. + ΔG_ele + ΔG_sol) is the binding energy calculated using the Generalized Born/Surface Area model (GB/SA) method using the last 25 ns of the trajectory (250 frames). The values in parenthesis are the standard mean error (SME) and reflect how precise the mean value is as an estimate of the true mean. The ligand-receptor binding energy (B.E.) is the binding energy calculated from docking. The following structural parameters are provided: distance between the C- and N-terminal (dCN), distance between ligand and C-terminal (dCL) and distance between ligand and N-terminal (dNL).

Pharmacophore modeling

The compounds were grouped according to the substitution position, namely G1-G4 (Figure 2). The first series includes the analogues of leucopleurotin, 2, 7, 8 and 9, which have a p-phenol group in their structure at G1. The analogues of dihydropleurotinic acid contain the compounds 3, 5 and 6, which have a carboxylic acid at G2. The other two sets are composed of the analogues of pleurogrisein, 4, 10 and 11, having a cyclopentane at G3 instead of an oxepane group, and the compounds 12, 13 and 14, which have a -SR group at G4 with a carboxylic acid at G2. Figure S2 (SI section) shows the optimized geometries with the molecular electrostatic potential (MEP) calculated for all 14 compounds.

The MD simulations were performed with the solvent (water) represented as a continuum model (SCRF approach). This is a satisfactory approximation for the present case, in which the binding site is a closed and deep cavity with low accessibility by the solvent. Moreover, in our previous study on pleurotin where explicit solvent was included, we did not observe the participation of solvent molecules on the ligand binding.¹⁶ The catalytic mechanism of TrxR involves the displacement of C-terminal arm towards the N-terminal inner moiety (see Figure 1). Therefore, the distances between C-terminal/N-terminal (dCN), C-terminal/ ligand (dCL) and N-terminal/ligand (dNL) are important parameters and were monitored along the MD trajectories. In addition, the average ligand-receptor binding energy was calculated from the last 25 ns of MD trajectory using GB/SA method (Table 1).³²

The mean square deviation (RMSD) in Table 1 shows that the complexes remain stable and with low structural fluctuations along the MD trajectory (ca. 1.5 Å), which is close to the free TrxR, 1.5 ± 0.3 Å. The arrangements of the ligands inside the binding site and their effects on the enzyme flexibility were analyzed. The binding mode expected for a potential inhibitor of TrxR is that in which the ligand remains stable between the N- and the C-terminal moieties, blocking their approach and forming a network of contacts with the enzyme (Figure 4). The distance dCN is 12.0 ± 0.9 Å for the free TrxR, and for the ligand-TrxR complexes this was longer, 13-22 Å (Table 1), except for compound 2, which presented dCN = 11 Å. In this case, the ligand moved to the inner part of the site, allowing the C-terminal arm approaches to N-terminal (minimum dCN = 7.83 Å, Figure 5). In general, there is a preference for the ligands remain closer to the C-terminal arm (dCL = 5-11 Å) than to the N-terminal (dNL = 8-12 Å).

Figure 5
Representation of ligand movement along the MD trajectory for the ligand 2.

Figure 6 shows the most suitable orientation of the ligands that is complementary to the catalytic site. The G1 moiety (p-quinone ring) is oriented towards the site mouth, whereas the G2 region (lactone ring) is pointed to the chain D (N-terminal) and G3 (oxepane ring) is oriented to the inner part of the site. This ligand-receptor arrangement favors the network of contacts needed to position the ligand between N- and C-terminal and, therefore, block the movement of C-terminal arm.

Figure 6
The most suitable conformation assumed by the ligands in the catalytic site (pleurotin 1 in this case). The surface of the catalytic site was constructed based on the residues part of the network of contacts. The red color represents the O atoms, the yellow the S atoms, blue the N atoms and the gray are the C atoms.

Regarding the stability of the complexes, the ligand-receptor binding energy was estimated from the MD trajectory using the GB/SA method. Table 1 shows the ΔG_aq, and in Table S1 (SI section), the main individual contributions are provided. The pleurotin 1 was taken as reference, ΔG_aq = -166 ± 21 kJ mol^-1, with more than 70% from van der Waals interactions. From Table 1, complex 2 was slightly more stable than pleurotin (–169 ± 14 kJ mol^-1) and complex 13 was the most stable among the 14 compounds evaluated, with GB/SA free energy of –211 ± 17 kJ mol^-1. The additional stability of 13 comes from the long and flexible -SR side chain in the p-quinone moiety, which allows more contacts with the residues at the catalytic site (Figure 7).

Figure 7
Representation of the network of contacts formed by compound 13 with TrxR. The main hot-spots are highlighted.

Table 2 includes the hot-spots, which are residues contributing with more than 4 kJ mol^-1 for the complex stability. The individual energy contributions are also shown in Table 2 and provide a quantitative picture of the weight of each residue to the complex stability. The complex 13 presents more hot-spots than the other ligands (12 considering both chains), which contribute with –79.1 kJ mol^-1 (37%) to the overall complex stability (see also Figure 7). Moreover, as predicted for pleurotin 1, the main contribution to the binding energy comes from van der Waals interactions (see Table S1), characteristic of hydrophobic ligands; however, ligand 13 forms H-bond with several polar residues (see Figure 7) that enhance the electrostatic contribution to 41%, compared to 27% for pleurotin. Using the last columns of Table 1 we note that ligand 13 has low score in docking, but high stability in MD simulation. This demonstrates the primary role of dynamics for flexible ligands.

Table 2
Hot-spots residues identified for the pleurotin derivatives in the binding site of TrxR

Hot-spots residues in chain C (ligand-residue binding energy / (kJ mol-1))
1	2	3	4	5	6	7	8	9	10	11	12	13	14
Thr471 (-6.95)	Phe397 (-9.33)	His463 (-4.64)	Leu400 (-5.03)	His463 (-6.87)	Thr472 (-9.16)	Phe397 (-7.17)	Ile483 (-5.78)	Leu484 (-6.90)	Phe397 (-5.99)	Leu400 (-4.44)	Val469 (-4.18)	Thr472 (-9.71)	Gly487 (-7.34)
Gln485 (-7.51)	Trp398 (-6.49)	Ile483 (-4.73)	His463 (-4.76)	Gln485 (-7.57)	Gln485 (-11.3)	Val469 (-9.70)	Gln485 (-6.01)	Gln485 (-13.7)	Pro399 (-4.54)	His463 (-9.30)	Thr472 (-7.10)	Leu473 (-6.71)	Cys488 (-8.64)
Cys489 (-7.10)	Leu400 (-6.11)	Gly487 (-6.23)	Val469 (-4.19)	Cys488 (-4.95)	Cys488 (-6.97)	Cys488 (-4.27)			His463 (-7.05)	Glu468 (-5.11)	Ile483 (-10.2)	Ser474 (-6.00)
	Ala467 (-5.59)				Cys489 (-5.02)				Val469 (5.21)	Thr472 (-4.28)		Val475 (-5.08)
	Glu468 (-8.06)								Thr472 (-9.85)	Ala486 (-6.21)		Gln485 (-10.2)
	Val469 (-6.94)								Gln485 (-5.56)	Cys489 (-5.55)		Ala486 (-4.85)
	Leu484 (-5.28)								Ala486 (-5.59)			Cys489 (-9.74)
	Gln485 (-5.66)								Gly485 (-8.64)
									Cys488 (-6.74)
Hot-spots residues in chain D (ligand-residue binding energy / (kJ mol-1))
1	2	3	4	5	6	7	8	9	10	11	12	13	14
Ala510 (-4.63)		Leu509 (-5.72)	Lys513 (-6.35)	Lys513 (-5.56)	Lys513 (-6.49)	Val544 (-5.18)	Leu509 (-4.21)	Leu509 (-5.10)	Lys552 (-5.24)	Ile549 (-9.13)	Lys513 (-4.86)	Leu509 (-4.49)	Val544 (-4.56)
Ile831 (-8.14)		Ala510 (-4.79)	Ile831 (-5.73)	Val544 (-6.07)	Ile831 (-5.55)	Ile831 (-4.67)	Ala510 (-4.63)	Leu544 (-6.40)	Ile831 (-4.49)	Ile831 (-4.33)	Ile831 (-5.05)	Ala510 (-5.64)	Tyr600 (-11.1)
		Lys513 (-5.44)	Arg835 (-6.03)	Ile831 (-7.70)	Arg835 (-6.62)		Lys513 (-1.60)	Ilr831 (-7.05)				Lys513 (-5.56)	Ile831 (-5.09)
		Tyr600 (-6.83)					Ile831 (-6.69)					Val544 (-5.13)
		Ile831 (-5.05)										Ile831 (-6.07)
(-34.3)	(-53.5)	(-43.4)	(-32.1)	(-38.7)	(-51.1)	(-30.9)	(-28.9)	(-39.1)	(-68.9)	(-48.3)	(-31.3)	(-79.1)	(-36.7)

• Residues with a binding energy contribution lower than -4 kJ mol^-1. In bold are the residues with the lowest energy among the hot-spots for each complex. The individual energy contributions are provided in parenthesis, and the total contribution is included in the last line of the table.

The data in Table 1 were used to analyze the role of substituents in the pleurotin backbone for the ligand-receptor interaction energy. Within the groups of molecules, it was observed that the ligands substituted at G2 (3, 5, 6) showed similar ΔG_aq, being ca. 10 kJ mol^-1 greater than for pleurotin 1 (see Table 1). At first glance, it suggests that the opening of lactone ring decreases the ligand-receptor interactions. To confirm this finding, pleurotin 1 can be compared with compound 3, and the compound 2 with 7, which pairs differ in the opening of the lactone ring at G2. As shown in Table 1, for the first pair (1/3) ΔG_aq increases by 7 kJ mol^-1 and for the second pair (2/7) the increasing was even larger, ca. 36 kJ mol^-1 due the opening of the lactone ring.

When comparing the lactone and the acid groups at G2, we observed that the acid interacts with Gly540 by H-bond. This interaction contributes to reorientate the ligand in the binding site, making it less favorable due the loss of several van der Waals contacts (see Table S1). Thus, it is possible to conclude that the lactone group at G2 favors the interaction with TrxR compared to the acid form.

When it comes to G3 region (ligands 4, 10, 11), it was observed that the molecules have different binding modes to TrxR, with a gap in the ΔG_aq of ca. 30 kJ mol^-1 between the most (10) and the least stable derivative (11). In general, it was observed that the opening of oxepane increases the degrees of freedom of the ligand, disrupting the favorable network of contacts that anchors the ligand (see Figure 4 to refer to the network of contacts); therefore, the drug-receptor complex was less stable than the parent compound pleurotin. Nonetheless, the ligand 10 presents ΔG_aq (= -163 kJ mol^-1) close to pleurotin, which is due to the conformation assumed during the MD simulation. Unlike the other molecules in G3, the compound 10 moved in the binding site and enhanced the interactions with the enzyme (Figure S3, SI section).

The molecular regions labeled as G1 and G4 present variations on p-quinone ring. The molecules 2, 7-9, have a reduced ring and 12-14 have a -SR substituent at the ortho position. The difference between the interactions with the TrxR of a p-phenol and p-quinone derivatives was analyzed, in conjunction to the role of the -SR substituent. Regarding the G1 position, except compound 2, the other analogues (7-9) with p-phenol group were less stable than the pleurotin 1, allowing us to infer that the oxidized form p-quinone should be part of the pharmacophore model. For compound 2, as represented in Figure 5, the ligand moved along of C-terminal arm and remains at the inner part of the binding site, stabilizing the complex. In spite of stronger interaction, the position of the ligand should not avoid the N- and C-terminal approach. To evaluate the p-quinone role for the pharmacophore model, compounds 3 and 7 can be compared. According to Table 1, the ΔG_aq indicates that complex 3 was more stable than 7, keeping the suitable conformation for binding TrxR. This arrangement has the hydrophilic region (labeled as G1 and G4 in Figure 2) facing the mouth opening of the binding site and the hydrophobic portion (G3) exposed directed towards the inner part of the cavity. Moreover, from Table 2, it is possible to observe that ligand 3 presents more hot-spots (eight) than ligand 7 (five) suggesting more effective contacts of 3 with TrxR. In addition to the number of host spots, their energy contributions were higher for binding of ligand 3 (-43.4 kJ mol^-1) than ligand 7 (-30.9 kJ mol^-1) (last line of Table 2).

The compounds 12, 13 and 14 presented ΔG_aq = –155 ± 21, -211 ± 17 and -155 ± 13 kJ mol^-1, respectively. The complexes 12 and 14 presented stability similar to complex 3 (ΔG_aq = -159 ± 15 kJ mol^-1) and smaller than pleurotin, which suggest that the small -SR substituent at the G4 did not increase the interactions with the catalytic site for those analogues. Conversely, for complex 13, which has a long -SR side chain, we predicted the lowest binding energy among the 14 ligands evaluated, -211 ± 17 kJ mol^-1. This molecule is similar to compound 3, with a long and flexible -SR side chain at G4 (p-quinone) that improved the contacts with the target, stabilizing the complex by ca. 60 kJ mol^-1 (see Figure 7).

The previous analysis allowed us to propose a minimum pharmacophore based on the pleurotin structure. The potential TrxR inhibitors should have the lactone group at G2, p-quinone at G1/G4 and an oxepane at G3 (see Figure 2 to refer to these positions), which corresponds to the basic pleurotin backbone (1) (Figure 8). The pharmacophore model represented in Figure 8 was constructed based on the main favorable interactions with TrxR observed for the ligands tested. When it comes to the residues that anchor the ligand, the His463, Val469, Thr472, Gln485, Cys488, Lys513 and Ile831 should be emphasized, once they were found on five or more ligands as hot-spots (see Table 2). His463 has an imidazole group that interacts through π-alkyl contact with the non-polar group at G3 (oxepane ring), thus it assists on the adjustment of the ligand in the site.Gln485 is in the middle of C-terminal arm and contributes to the alignment of the arm with the ligand. The charged Lys513 is located at the mouth opening of the binding site and acts like a “barrier” to ligand moves away the site. Moreover, it forms a cation-π interaction with the p-quinone group at G1. The Tyr600 acts as H-bond donor to the carbonyl groups at G2 and G1, although with lower frequency than the other hot-spots. However, Tyr600 and His463 have a primary role to the catalytic mechanism of TrxR,¹⁷ whose interactions with ligand may prevent the enzyme action. Lastly, the Ile831 was highlighted as a hot-spot for all complex, except for complex 2. This residue interacts with the ligands by van der Waals forces through the alkyl side chain of the residue with the G1 and G3 portions.

Figure 8. Pleurotin-based pharmacophore for inhibitors of TrxR. The naturally occurring pleurotin backbone must be retaining and the side groups R¹, R² and R³ are used to enhance the TrxR affinity.

Figure 8
Pleurotin-based pharmacophore for inhibitors of TrxR. The naturally occurring pleurotin backbone must be retaining and the side groups R1, R2 and R3 are used to enhance the TrxR affinity.

Design of new pleurotin analogues

The pharmacophore model represented in Figure 8 was used to design 20 new analogues with substitutions at R¹, R² and R³ (see Figure 3). The new structures were firstly subjected to docking simulations with the very same protocol used for the set of tests. The B.E. are given in Table S2 and Figure S1 (SI section). Nine ligands presented B.E. ≤ -40 kJ mol^-1, being the best representatives ligands with R³ = -CH₂NH₂, namely, 22, 32, 33 and 34 (Figure 9).

Figure 9
Interactions between TrxR and the compounds 22, 32, 33 and 34 obtained from molecular docking. The highlighted residues improve the interactions.

Taking the pleurotin as a reference, it was observed that the non-polar linear groups (R³ = –Et (15), –Pr (21)) slightly increase the complex stability allowing the ligand to access the internal hydrophobic part of the binding site. Conversely, the non-polar and bulky groups (R³ = –i-Pr (16) –t-butyl (17)) decrease the complex stability due to the high steric hindrance. These bulky ligands reduce the interactions with the polar amino acids His463 and Glu468 at the catalytic site. When hydroxide groups are included at R³ (18, 19, 20), the stability decreases significantly. This effect is due to the increased polarity on the hydrophobic part of the pharmacophore (G3, oxepane ring). The analogues with R³ = –CH₂NH₂ (22, 32, 33 and 34) favor H-bond with Glu468 and Cys488 as seen in Figure 9. When it comes to R¹ and R² positions, the groups –Me, –OH, –Cl, –NH₂ were tested, keeping R³= –CH₃. In general, it was noted that small substituents at R² favor the interaction, whereas substitution in R¹ results in a decrease of the negative charge near the carbonyl moieties that destabilizes the complex by preventing interactions of the carbonyl of the p-quinine ring with the binding site.

By comparing the MEP for the pairs 26/30, 24/28 and 27/31 (Figure S4, SI section), it is observed that substitution in R¹ enhances the MEP (more positive) in the p-quinone ring, which disfavors the ligand-receptor interactions. Besides, the –CH₃ group in R² position presented the best score, increasing the interactions with the C-terminal and fixing the ligand in a favorable position within the network of contacts. The previous conclusions based on the ligands 15-31, suggest the following structures as the best representatives: 32 (R¹ = –H, R² = –CH₃ and R³ = –CH₂NH₂), 33 (R¹ = –SCH₂CH(NHCOCH₃)COOH, R² = –CH₃ and R³ = –CH₂NH₂) and 34 (R¹ = –SCH₂CH(NHCOCH₃) COOH, R² = –H and R³ = –CH₂NH₂). The docking scores were B.E. = -44.14, -43.10 and -42.51 kJ mol^-1, respectively, being 6-10% more stable than pleurotin.

To improve and validate the results, MD simulation was carried out for the new analogues 32, 33 and 34. The GB/SA ΔG_aq was -178 ± 15, -236 ± 21 and -213 ± 21 kJ mol^-1, respectively, found more stable than the reference pleurotin 1 (ΔG_aq = -166 ± 21 kJ mol^-1) and the analogue 13 (-211 ± 17 kJ mol^-1) from the original set. It is worth noting that compounds 13 and 33 (most stable) differ by the substituent at R² and R³ (-CH₃ and -CH₂NH₂, respectively) and by the lactone ring at G2. The dCN was 17 ± 1 (32), 13.9 ± 0.7 (33) and 19 ± 1 Å (34), indicating that these compounds were able to block the approach of N- and C-terminal, which is necessary to inhibit the TrxR function. Figure 10 shows the equilibrium conformations assumed for complexes 32, 33 and 34 where the main residues with frequency of contacts > 40% are highlighted. For the best ligand 33, the Cys488 (chain C), Lys513 (chain D) and Tyr600 (chain D) play a primary role for complex stability, interacting with the three anchor points identified in the pharmacophore, i.e., R³ = -CH₂NH₂, R¹ = -SR (long and flexible side chain) and the lactone ring at G2 (see also Figure 8).

Figure 10
Structures of the TrxR complexes with the ligands 32, 33 and 34. The highlighted residues were selected from MD trajectory based on the frequency of contacts > 40%. The B.E. is the binding energy calculated from docking, and the ΔG_aq were calculated with GB/SA method from the last 25 ns of MD trajectory.

In summary, the pleurotin-based pharmacophore model proposed here lead to a potential representative as TrxR inhibitor, labeled as compound 33. This molecule was not synthesized yet but is structurally similar to the natural analogue 13, which could be used as a lead compound for a semisynthetic route of pleurotin derivatives. Moreover, this work will certainly motivate the experimentalists to engage in the synthesis of compound 33 and other similar analogues. Finally, it is worth mentioning that pleurotin and its derivatives have biological potential that goes beyond the inhibition of TrxR, therefore, the design and synthesis of new analogues are of continuing interest, and the present analysis and results may help to advance in this way.

Conclusions

A pharmacophore for inhibitors of TrxR was proposed based on natural analogues of pleurotin, which is known to inhibit the TrxR enzyme with IC₅₀ = 0.17 μM. Docking and molecular dynamics simulations were used to monitor the ligand-TrxR interactions and to assess the main structural features necessary for enzyme binding. We found out that the pleurotin backbone must be preserved, but side substituents at key positions can be included to enhance the interactions. The ligand-TrxR binding free energy (∆G_aq / kJ mol^-1) was -166 and decreased to -211 with the -SR condensation on the p-benzoquinone moiety (ligand 13). Among the 12 TrxR residues assigned as hot-spots for ligand 13, the Thr473, Gln485 and Cys489, in the C-terminal domain, contributed mostly to the binding energy, due to the interaction with the -SR side group. The pharmacophore model proposed for TrxR inhibitors was used to propose new plurotion derivatives, and compound 33, which is analogue to 13, differing by the lactone ring at G2, –CH₃ at R² and –CH₂NH₂ at R³, formed a very stable complex with TrxR (∆G_aq = –236 kJ mol^-1). In addition to the network of contacts found for ligand 13, the ligand 33 interacted with the Gly487 through H-bond with the –CH₂NH₂ substituent, which favored others favorable contacts with Lys513 and Tyr600. This ligand was not synthetized yet, but might serve as lead compound for further experimental studies. Moreover, the pleurotin-based pharmacophore presented here can also be used to design pleurotin derivatives with biological potential beyond TrxR inhibitors.

Supplementary Information

Supplementary data are available free of charge at http://jbcs.sbq.org.br as PDF file.

Acknowledgments

The authors would like to express gratitude to the financial support from the Brazilian Agencies CNPq (307018/2021-0), FAPEMIG and CAPES that provided continuous support to our laboratories.

References

• ¹
  International Agency for Research on Cancer (IARC), Key Cancer Data and Key Figures on IARC: 2020-2021, https://www.iarc.who.int/biennial-report-2020-2021web, accessed in January 2024.
  » https://www.iarc.who.int/biennial-report-2020-2021web,
• ²
  Lee, Y. T.; Tan, Y. J.; Oon, C. E.; Eur. J. Pharmacol. 2018, 834, 188. [Crossref]
  » Crossref
• ³
  Karunanithi, S.; Liu, R.; Hou, Y.; Gonzalez, G.; Oldford, N.; Roe, A. J.; Idipilly, N.; Gupta, K.; Amara, C. S.; Putluri, S.; Lee, G. K.; Valentin-Goyco, J.; Stetson, L.; Moreton, S. A.; Putluri, V. ; Kavuri, S. M.; Saunthararajah, Y. ; de Lima, M.; Tochtrop, G. P.; Putluri, N.; Wald, D. N.; Oncogene 2021, 40, 5236. [Crossref]
  » Crossref
• ⁴
  Bian, M.; Fan, R.; Zhao, S.; Liu, W.; J. Med. Chem. 2019, 62, 7309. [Crossref]
  » Crossref
• ⁵
  Zhang, J.; Li, X.; Han, X.; Liu, R.; Fang, J.; Trends Pharmacol. Sci. 2017, 38, 794. [Crossref]
  » Crossref
• ⁶
  Zhang, J.; Zhang, B.; Li, X.; Han, X.; Liu, R.; Fang, J.; Med. Res. Rev. 2019, 39, 5. [Crossref]
  » Crossref
• ⁷
  Lu, J.; Holmgren, A.; Free Radicals Biol. Med. 2014, 66, 75. [Crossref]
  » Crossref
• ⁸
  Zhang, X.; Selvaraju, K.; Saei, A. A.; D’Arcy, P.; Zubarev, R. A.; Arnér, E. S. J.; Linder, S.; Biochimie 2019, 162, 46. [Crossref]
  » Crossref
• ⁹
  Hashemy, S. I.; Ungerstedt, J. S.; Avval, F. Z.; Holmgren, A.; J. Biol. Chem. 2006, 281, 10691. [Crossref]
  » Crossref
• ¹⁰
  Harper, M. T.; Platelets 2019, 30, 98. [Crossref]
  » Crossref
• ¹¹
  Qiu, X.; Liu, Z.; Shao, W.-Y.; Liu, X.; Jing, D.-P.; Yu, Y.-J.; An, L.-K.; Huang, S.-L.; Bu, X.-Z.; Huang, Z.-S.; Gu, L.-Q.; Bioorg. Med. Chem. 2008, 16, 8035. [Crossref]
  » Crossref
• ¹²
  Sweeney, M.; Coyle, R.; Kavanagh, P.; Berezin, A. A.; Lo Re, D.; Zissimou, G. A.; Koutentis, P. A.; Carty, M. P.; Aldabbagh, F.; Bioorg. Med. Chem. 2016, 24, 3565. [Crossref]
  » Crossref
• ¹³
  Powis, G.; Wipf, P.; Lynch, S. M.; Birmingham, A.; Kirkpatrick, D. L.; Mol. Cancer Ther. 2006, 5, 630. [Crossref]
  » Crossref
• ¹⁴
  Hoskin, J. F.; Sorensen, E. J.; J. Am. Chem. Soc. 2022, 144, 14042. [Crossref]
  » Crossref
• ¹⁵
  Wipf, P.; Hopkins, T. D.; Jung, J.-K.; Rodriguez, S.; Birmingham, A.; Southwick, E. C.; Lazo, J. S.; Powis, G.; Bioorg. Med. Chem. Lett. 2001, 11, 2637. [Crossref]
  » Crossref
• ¹⁶
  Quintanilha, D. B.; Dos Santos, H. F.; J. Biomol. Struct. Dyn. 2023, 41, 5646. [Crossref]
  » Crossref
• ¹⁷
  Fritz-Wolf, K.; Urig, S.; Becker, K.; J. Mol. Biol. 2007, 370, 116. [Crossref]
  » Crossref
• ¹⁸
  Er-rajy, M.; Fadili, M. E.; Mujwar, S.; Lenda, F. Z.; Zarougui, S.; Elhallaoui, M.; Struct. Chem. 2023, 34, 1527. [Crossref]
  » Crossref
• ¹⁹
  Sandargo, B.; Thongbai, B.; Praditya, D.; Steinmann, E.; Stadler, M.; Surup, F.; Molecules 2018, 23, 2697. [Crossref]
  » Crossref
• ²⁰
  Sandargo, B.; Thongbai, B.; Stadler, M.; Surup, F.; J. Nat. Prod. 2018, 81, 286. [Crossref]
  » Crossref
• ²¹
  Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P.; J. Am. Chem. Soc. 1985, 107, 3902. [Crossref]
  » Crossref
• ²²
  Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V. ; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y. ; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V. ; Cioslowski, J.; Fox, D. J.; Gaussian 09, Revision D.01; Gaussian, Inc.: Wallingford, CT, 2004.
• ²³
  Protein Data Bank (PDB), https://www.rcsb.org/structure/2J3N, accessed in January 2024.
  » https://www.rcsb.org/structure/2J3N,
• ²⁴
  Humphrey, W.; Dalke, A.; Schulten, K.; J. Mol. Graphics 1996, 14, 33. [Crossref]
  » Crossref
• ²⁵
  PyMOL, version 1.8; Schrödinger LLC, New York, U.S.A., 2015.
• ²⁶
  Case, D. A.; Betz, R. M.; Cerutti, D. S.; Cheatham III, T. E.; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.; Homeyer, N.; Izadi, S.; Janowski, P.; Kaus, J.; Kovalenko, A.; Lee, T. S.; LeGrand, S.; Li, P.; Lin, C.; Luchko, T.; Luo, R.; Madej, B.; Mermelstein, D.; Merz, K. M.; Monard, G.; Nguyen, H.; Nguyen, H. T.; Omelyan, I.; Onufriev, A.; Roe, D. R.; Roitberg, A.; Sagui, C.; Simmerling, C. L.; Botello-Smith, W. M.; Swails, J.; Walker, R. C.; Wang, J.; Wolf, R. M.; Wu, X.; Xiao, L.; Kollman, P. A.; AMBER 2016, University of California, San Francisco, 2016.
• ²⁷
  Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J.; J. Comput. Chem. 2009, 30, 2785. [Crossref]
  » Crossref
• ²⁸
  Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C.; J. Chem. Theory Comput. 2015, 11, 3696. [Crossref]
  » Crossref
• ²⁹
  Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A.; J. Comput. Chem. 2004, 25, 1157. [Crossref]
  » Crossref
• ³⁰
  Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A.; J. Mol. Graphics Model. 2006, 25, 247. [Crossref]
  » Crossref
• ³¹
  Uberuaga, B. P.; Anghel, M.; Voter, A. F.; J. Chem. Phys. 2004, 120, 6363. [Crossref]
  » Crossref
• ³²
  Miller, B. R. I.; McGee Jr., T. D.; Swails, J. M.; Homeyer, N.; Gohlke, H.; Roitberg, A. E.; J. Chem. Theory Comput. 2012, 8, 3314. [Crossref]
  » Crossref

Edited by

Editor handled this article: Paula Homem-de-Mello (Associate)

Publication Dates

• Publication in this collection
  12 Feb 2024
• Date of issue
  2024

History

• Received
  31 Oct 2023
• Accepted
  31 Jan 2023