Unveiling the Possible Oryzalin-Binding Site in the α-Tubulin of Toxoplasma gondii.

Dinitroaniline derivatives have been widely used as herbicidal agents to control weeds and grass. Previous studies demonstrated that these compounds also exhibit good antiparasitic activity against some protozoan parasites. Oryzalin (ORY), a representative dinitroaniline derivative, exerts its antiprotozoal activity against Toxoplasma gondii by inhibiting the microtubule polymerization process. Moreover, the identification of ORY-resistant T. gondii lines obtained by chemical mutagenesis confirmed that this compound binds selectively to α-tubulin. Based on experimental information reported so far and a multiple sequence analysis carried out in this work, we propose that the pironetin (PIR) site is the potential ORY-binding site. Therefore, we employed state-of-the-art computational approaches to characterize the interaction profile of ORY at the proposed site in the α-tubulin of T. gondii. An exhaustive search for other possible binding sites was performed using the Wrap "N" Shake method, which showed that ORY exhibits highest stability and affinity for the PIR site. Moreover, our molecular dynamics simulations revealed that the dipropylamine substituent of ORY interacts with a hydrophobic pocket, while the sulfonamide group formed hydrogen bonds with water molecules at the site entrance. Overall, our results suggest that ORY binds to the PIR site on the α-tubulin of the protozoan parasite T. gondii. This information will be very useful for designing less toxic and more potent antiprotozoal agents.

Introduction

Oryzalin (ORY, Figure A) is a pre-emergent dinitroaniline herbicide widely used to control annual grasses and broadleaf weeds on fruit trees and crops.[1] Like trifluralin (TFL, Figure A) and other dinitroaniline derivatives, ORY exerts herbicidal activity by selectively disrupting plant microtubules.[2,3] Interestingly, this compound has also been shown to inhibit microtubule polymerization of protozoan parasites such as Toxoplasma gondii, Plasmodium falciparum, Leishmania mexicana spp., and Trypanosoma cruzi.[4−8] The low effect on mammalian tubulin polymerization and low inhibitory concentrations required to kill the parasites have pinpointed ORY as a lead compound for the development of new antiprotozoal compounds.[5,7,9−11]

Figure 1

(A) Chemical structures of the dinitroaniline ORY (left) and TFL (right). (B) Depiction of dinitroaniline resistance mutations (purple spheres) reported in the α-tubulin of protozoan parasites and plants. (C) Distribution of the 30 mutations reported in TgAT colored by the OR concentrations shown at the bottom of the figure. (D) Schematic representation of the L-, M-, and N-loops in the α-tubulin monomer and residues comprising the GTP binding site (GTPbs).

A total of 30 mutations associated with ORY resistance (OR) in 23 different amino acids have been reported in T. gondii α-tubulin (TgAT, Figure B).[12−17] Dinitroaniline resistance mutations in plants and T. gondii suggest that the dinitroaniline binding site could be located in the α-tubulin.[14,16−18] Morrisette and co-workers[13] determined that most of these mutations confer low (OR = 0.5–2.0 μM) to moderate (OR = 2.0–25 μM) levels of OR when compared to the wild-type parasite (IC50 = 0.1 μM).[4] Nevertheless, TgAT with L136F, T239I, R243S, or V252L mutation was resistant to higher concentrations of ORY (OR > 35 μM[13]) (Figure C).

Based on the OR mutation data and computational results, several research groups have proposed different ORY binding sites (ORYbs) in α-tubulin of plants and protozoan parasites. Some of these groups have suggested that dinitroaniline herbicides could bind directly behind the α-tubulin N-loop (H1–S2 loop, amino acids 29–64) or close to the L-loop (S9–S10 loop, amino acids 356–373).[14,19−21] Meanwhile, other studies have proposed that ORYbs might be located on top of the N-loop near to the β-tubulin interface (see Figure D).[22−24] Unfortunately, these proposed binding sites are found in poorly accessible regions, away from amino acids with known OR mutations or comprise several highly conserved residues between ORY-susceptible and non-susceptible organisms.

Recently, Chen et al.[25] have found that V202F α-tubulin mutation conferred dinitroaniline resistance to Lolium rigidum seedlings, as previously reported by Hashim et al.[26] with the perennial grass Alopecurus aequalis. Qin Yu group also carried out a multiple sequence alignment of α-tubulins in which they discovered that non-susceptible organisms have a Phe residue at position 202. In contrast, highly susceptible organisms (i.e., protists and plants) have Val or Ile at this position. This remarkable observation revealed that this residue could play a critical role in the interaction of dinitroaniline derivatives with α-tubulin, shedding new insights into the location of ORYbs.

Several microtubule-stabilizing agents and polymerization inhibitors have been co-crystallized with the tubulin.[27,28] Until the writing of this paper, unsaturated δ-lactone pironetin (PIR) was the only tubulin polymerization inhibitor that binds to the α-tubulin core (Figure A).[29] Interestingly, most of the mutations that have been reported so far are located within this α-tubulin region.[30] Therefore, we propose that ORY binds to the PIR site in α-tubulin. To prove this, we collected the most recent experimental data regarding ORY tubulin polymerization inhibitory activity in protozoan parasites and plants. Subsequently, ORY was docked into TgAT using molecular docking and Wrap “N” Shake (WnS) protocols. The best ranked TgAT–ORY complex was further subjected to three independent 200 ns molecular dynamics (MD) simulations. Finally, the five most representative complexes were extracted to compute ORY binding free energies using alchemical calculations.

Figure 2

(A) Depiction of the PIRbs in the Sus scrofa α-tubulin (SsAT) (PDB ID: 5fnv). (B) PIR chemical structure and representation of residues involved in PIR interaction. (C) Final ORY pose in the PIRbs of TgAT from the flexible receptor molecular docking studies. Hydrogen bonds are represented as red dotted lines. The flexible residues are shown in black.

Results and Discussion

OR Mutation Residues Are Found at the Pironetin Site

It is well known that dinitroaniline derivatives exhibit herbicidal and antiprotozoal activity selectively binding to α-tubulin.[25] The tubulin polymerization inhibitory activity of these compounds has been observed in plants such as A. aequalis (water foxtail),[26]L. rigidum (annual grass),[25]Eleusine indica (goosegrass),[31] and Setaria viridis (green foxtail).[32] Similarly, treatment with dinitroaniline derivatives has shown microtubule disruption activity in protozoa with low half-maximal inhibitory concentrations (IC50). T. gondii (IC50 = 0.1 μM),[4]P. falciparum (IC50 = 3.9 μM),[33] and L. mexicana (IC50 < 20 μM)[7,11] are some of the most relevant ORY-susceptible protozoan parasites. A previous study suggested that the protozoan Giardia intestinalis was also susceptible to ORY treatment.[11] However, higher concentrations of this dinitroaniline are required to kill this parasite (IC50 > 50 μM). Something similar occurs with the nematode model Caenorhabditis elegans and in mammalian cell lines, where cytotoxic activities have been observed at IC50 values greater than 25 and 50 μM, respectively.[4,5,9,34,35] In agreement with the conclusions made by Terra et al.[11] and Hirst et al.,[33]G. intestinalis, C. elegans, and mammals could be considered ORY-non-susceptible organisms.

To understand the difference in the susceptibility to ORY in these organisms, we performed a multiple sequence alignment of the full-length α-tubulin amino acid sequences (Figure S1 of the Supporting Information). For this analysis, we included α-tubulin sequences of disease-causing nematodes Trichinella spiralis (trichinosis), Ascaris suum (ascariasis), and Haemonchus contortus (haemonchosis). Some of the most important changes between the sequences of the susceptible and non-susceptible organisms were found at positions 200, 202, 238, 252, 268, and 353. Variations at position 202 turn out to be one of the most interesting changes among α-tubulin sequences because they are consistent with observations made by the Qin Yu group.[25] In that study, the authors mentioned that replacing a small residue (i.e., V or I) at position 202 with a phenylalanine residue leads to dinitroaniline resistance. A recent structural analysis carried out by Gaillard et al.[36] with the protozoan Tetrahymena thermophila α/β-tubulin heterodimer has revealed that V202 is located at the PIR binding site (PIRbs). Despite this finding, little importance was given to characterize PIRbs because the focus has been centered in the identification of new tubulin polymerization inhibitors that bind to the colchicine/benzimidazole site on β-tubulin. It is worth mentioning that PIRbs on α-tubulin is in the same region as the benzimidazole site (BZbs) on β-tubulin. Contrary to ORY, benzimidazole derivatives do not bind to β-tubulin of T. gondii, P. falciparum, and L. mexicana but do interact with β-tubulin of G. intestinalis, nematodes, and mammals.[37]

In addition, our multiple sequence analysis suggests that most OR mutations in T. gondii constitute the PIRbs. A closer analysis showed that three of the four high resistance mutations (OR > 35 μM) are located at the PIRbs (L136F, T239I, and V252L) as well as five of the seven mutations with known moderate OR (V4L, F52L/Y, S165T/A/P, M268T, and I378M).[4,13] Overall, these observations support the hypothesis that ORY could occupy the same site as PIR on α-tubulin.

ORY Interacts with the Hydrophobic Pocket of PIRbs

PIR is a potent anticancer natural product that covalently binds to residue C316 of the α-tubulin subunit to arrest the cell cycle progression (Figure B).[30] This molecule has an α,β-unsaturated lactone responsible for covalent binding at the entrance of the site and an aliphatic chain that interacts with a deeper hydrophobic pocket.[29] To determine the binding mode of ORY in PIRbs of the TgAT model (Figure S2 of the Supporting Information), we employed different molecular docking approaches. Our rigid receptor molecular docking study showed that the dipropylamine chains of ORY are oriented toward the hydrophobic pocket of the PIRbs, while the sulfonamide substituent shares the same site as the PIR lactone ring. This result suggests that ORY binds to the proposed binding site and displays an interaction profile similar to PIR. To relax the interaction of the binding site residues with ORY, we performed several flexible molecular docking assays by testing different combinations of mobile side chains (Table S1 of the Supporting Information). Our study found that the mobility of L238, Q256, M268, C316, and M318 residues was key to improve the intermolecular interactions and docking score (Figure C). The final ORY binding mode displayed a hydrogen bond interaction between the nitro groups of the molecule and the main chain of L242 and the side chain of Q256. Furthermore, we found that the dipropylamine group of ORY was oriented in a more similar way to the pentenyl group of PIR in the hydrophobic pocket. These alkyl substituents formed hydrophobic interactions with non-polar side chains of L136, L167, L238, L242, and V252, while the aromatic ring of the dinitroaniline showed a single hydrophobic contact with the β carbon (Cβ) atom of F255. On the other hand, the sulfonamide remained at the pocket entrance without presenting any significant intermolecular interaction with the site residues. This hydrophilic substituent could instead be playing an essential role in the interaction of water molecules entering the site (vide infra).

To confirm the location of our proposed ORYbs and compare its binding affinity with other possible sites, we carried out an exhaustive search for other potential binding sites for this compound on α-tubulin using the WnS methodology (Figure ). We performed three independent cycles of the WnS protocol (WnS1, WnS2, and WnS3) to have more reliable results. First, we wrapped the entire TgAT structure with 198 to 226 copies of ORY molecules. The wrapped TgAT systems were subsequently subjected to five cycles of backbone-restrained MDs, one conventional (MDB), and four using a simulated annealing scheme (MDBSA). In the last MDBSA cycle, less than 24 copies of ORY remained bound to TgAT. ORY poses between replicates were compared and clustered to detect those binding sites in at least two of the three independent assays. In total, we found 15 potential ORY binding sites on α-tubulin (s01 to s15).

Figure 3

Schematic representation of the WnS workflow, number of ORY copies (NORY) bound to TgAT during the different shake cycles, and binding energies of the copies clustered at the 15 potential sites (s01 to s15) for the three WnS replicates (WnS1 to WnS3).

Nevertheless, low binding free energy values and pose consistency between the different replicas were obtained at only sites s01 to s03. Unfortunately, s02 and s03 sites are located far from residues associated with OR, precluding the possibility of being other potential ORY binding sites. Interestingly, the ORY copy in the PIRbs (s01) showed the lowest binding energy in all three replicates. In addition, the per-residue energy contribution analysis showed that s01 was the only binding site where several residues with known OR mutations play a key role in the interaction of TgAT with ORY (Figure S3 of the Supporting Information). This result confirms that this is the most likely ORYbs and that ORY has a higher affinity for this site than for any other pocket or cavity on TgAT, including the guanosine-5′-triphosphate site (GTPbs).

ORY Interacts with Key OR Mutation Residues

Our results demonstrated that ORY is highly likely to bind to the same site as PIR on TgAT. To understand the dynamic behavior of ORY at the PIRbs, we carried out three independent 200 ns MD simulations (MD1 to MD3) of the TgAT–ORY complex obtained from the flexible molecular docking study (Figure S4 of the Supporting Information). The GTP molecule and the Mg ion bounded to the GTPbs were also included during the simulations. The three MD replicas’ ligand root-mean-square deviation (RMSD) analysis showed that ORY suffered a conformational change of approximately 0.31 nm from the original pose during the simulated time (Figure A). Noteworthy, ORY did not leave the pocket in any MD simulations. Based on these results, the last 190 ns of each replicate was merged into a single trajectory (28,500 frames) to have a deeper and broader insight into the behavior of ORY in the PIRbs pocket. First, we compute the fraction of ORY contacts (QORY) with amino acids that constitute the ORYbs in TgAT (Figure B). This analysis showed that ORY interacted for more than 50% of the trajectory with residues L136 (64.6%), V202 (56.4%), T239 (69.5%), and V252 (89.8%). These residues stand out from the rest because they are four of the five amino acids whose modifications confer resistance to high concentrations of ORY in plants and protozoa.[16,23,25] Our analysis also showed that these key residues, along with other 11 amino acids that constitute the PIRbs, contact ORY through hydrophobic interactions (HI, Figure C). The interaction of these hydrophobic pocket residues occurs mainly with the ORY dipropylamine group. Furthermore, we found that ORY forms a stable π–π stacking (πS) interaction with F255 side chain during 42% of the MD trajectory (Figure D). This information, along with the per-residue energy contribution analysis (Figure E), demonstrates that ORY interacts with some of the most important residues whose mutations are experimentally known to generate resistance to treatment with this compound. On the other hand, hydrogen bonds between ORY and polar amino acids S237, S241, and N258 were observed for less than 25% of the whole trajectory, suggesting that these interactions are not essential for ligand stability (Figure F).

Figure 4

(A) RMSD of ORY heavy atoms after least square fit to TgAT backbone for the three independent 200 ns MD simulations (MD1 to MD3). (B) Fraction of the number of ORY contacts with the amino acids that constitute the proposed binding site in TgAT [f(QORY)]. The boxes at the top of the graph pinpoints the residues with reported OR mutations, colored as in Figure C. (C) Fraction of the number of hydrophobic interactions between ORY and the binding site residues [f(HI)]. (D) Simulated distributions of MD trajectory projected onto the distance between the center of geometry of ORY and F255 benzene groups (COGdist) and the angle formed between these aromatic rings (ω). TgAT-ORY complexes exhibiting a COGdist value lower than 4.5 Å and a ± 20° ω angle were counted as π–π stacking interactions (orange dots).[43] (E) Energy contribution of TgAT residues to ORY affinity calculated from the merged MD simulation. Residues with known OR mutations are colored in red. (F) Number of hydrogen bonds formed between ORY and the side chains of S237, S241, and N258 in the merged trajectory (28,500 frames). The numbers to the right of the graph show the hydrogen bond occupancy values for each interaction. (G) Number of contacts between water molecules and the (a) dipropylamine, (b) dinitrobenzene, and (c) sulfonamide segments of ORY (Qwater).

In addition to the hydrogen bonds listed above, we found that the ORY sulfonamide substituent forms many contacts with water molecules (Qwater) (Figure G). This behavior is not surprising due to its proximity to the binding site entrance and explains the orientation of the ligand on the site because these functional groups are known to have a high affinity for the aqueous solvent.[38] In contrast, fewer contacts between water molecules and nitro groups occurred during the simulations.

Identification of Key Structural Water at the ORYbs

During the MD analysis, we noticed the presence of a structural water molecule that bridges TgAT with one of the nitro groups (Figure A). This water molecule was found at a small hydrophilic cavity, forming hydrogen bonds with the main chain of S237, L238, S241, and L242. Surprisingly, this water molecule was also identified in a Bos taurus α-tubulin (BtAT) crystal structure complexed with PIR.[39] In our simulations, the structural water molecule interacts with L238, S241, and L242 residues during more than 60% of the trajectory (Figure B). Although different interaction patterns between the protein main chain and the ligand were evaluated, we found that the TgAT-ORY water bridge occurs in less than 20% of the simulated time. Despite this, this water molecule could play an essential role because something very similar has been seen in the interaction of ligands with the BZbs in the β-tubulin.[40] Previous studies have explored the importance of this structural water in BZbs, leading to the conclusion that this molecule stabilizes the H7 helix break and acts as a hydrogen bond donor group to interact with tubulin polymerization inhibitors.[40,41] Moreover, the presence of this water molecule in recent high-resolution crystallographic structures has become more relevant to determine the interaction profile of new destabilizing agents.[42]

Figure 5

(A) Depiction of the water site located at the ORY binding pocket in TgAT. The five most representative TgAT–ORY complexes of the merged trajectory were superimposed for the visualization of the water site. The upper right figure shows the interaction of the structural water molecule with ORY in TgAT, while the bottom right shows a similar solvent molecule in the B. taurus α-tubulin (BtAT) crystal structure complexed with PIR (chain C of PDB ID: 5LA6(39)). (B) Fraction of trajectory frames in which water bridges were formed between residues (AA–H2O) and/or with ORY (AA–H2O–ORY).

ORY Exhibits High Dynamic Stability and Affinity for PIRbs

As mentioned above, ORY did not leave the pocket on TgAT or show significant displacement from its original position. By examining the distribution of the center of geometry of the benzene nucleus (COGbz) and its angle with respect to the initial pose (θ), we found that ORY presented small displacements in the pocket through the trajectory with a subtle swinging movement (Figure A). RMSD-based clustering analysis, in combination with the two-dimensional distribution of COGbz and θ parameters, showed that there are five main conformational groups (clusters C01 to C05) that represent 92.7% of the trajectory (Figure B). Of these clusters, we noted that the conformations grouped in C03 are those with a very similar pose to the initial ORY binding mode. In the RMSD plot of the independent simulations (see Figure A), these conformers can be found in those frames where an RMSD ≤ 0.2 nm is reached. The most representative poses of each of these clusters (p01 to p05) were extracted from the trajectory to better understand the dynamic of ORY at the site (Figure C). From these structures, we could visualize the swinging movement of ORY and the hydrogen bonding interactions of its sulfonamide substituent with N258 (p01) and S241 (p04) and the predilection of the aliphatic chains of the dipropylamine group for the hydrophobic pocket. For each of the poses, the absolute binding free energy (ΔGbind) of the ligand was computed using the free energy perturbation (FEP) method (Figure D). The use of different initial positions of ORY allowed us to determine that the binding free energy of this compound is −9.82 ± 1.26 kcal/mol. As expected, very similar ΔGbind values were obtained between the different poses. The highest energy value was achieved for p01 and the lowest for p02, while the other three showed no significant changes.

Figure 6

(A) Schematic representation of the distribution of ORY benzene center of geometry (COGbz) during the trajectory colored by the angle between ring planes relative to the initial pose (θ). The figures at the bottom of the image show the parameters used for the analysis. (B) COGbz distribution in X and Y planes (2D) colored by θ (top) and the cluster number (middle). The bottom plot shows the fraction of frames that constitute the five clusters with larger conformational population (C01 to C05). (C) Depiction of the ORY binding mode from the flexible molecular docking (black) and the most representative poses (p01 to p05) extracted from the five clusters. (D) Non-physical thermodynamic cycle used in the FEP method for the absolute binding free energy calculation (ΔGbind) of ORY in the TgAT complex (right). Equation used to calculate the ΔGbind (upper right) for each of the ORY poses (lower right).

Our MD results confirm the stability and binding affinity of ORY for the PIRbs of TgAT. Unlike molecular docking approaches, the use of explicit water molecules during the simulations allowed us to understand the role of the sulfonamide substituent in the protein–ligand interaction. The nitro groups, along with the sulfonamide, seem to block water molecules’ internalization into the hydrophobic pocket. The influence of the dipropylamine group in the interaction with the hydrophobic pocket residues was also corroborated with these simulations. Surprisingly, some of these findings have also been observed at the interaction between PIR and tubulin, despite its significant structural difference from ORY.[39]

Conclusions

In this work, several computational approaches were used to identify and characterize the possible ORY-binding site on TgAT. With a multiple sequence alignment analysis, we found that the most critical differences in α-tubulin from susceptible and non-susceptible organisms were located at the pironetin site. Moreover, three of the four high ORY resistance mutations were located at the PIRbs as well as V202F (reported in plants) and five of the seven mutations with known moderate OR. Our rigid and flexible molecular docking studies showed that the aliphatic chains of the dipropylamine interact with a hydrophobic pocket, while the sulfonamide substituent points toward the entrance of the binding site. Despite the exhaustive search for other possible binding sites, ORY exhibited a greater affinity for the PIRbs. MD simulations of the TgAT-ORY complex showed that ORY has high stability in the complex and affinity for the binding site, presenting small conformational changes guided by the hydrophobic contacts of the dipropylamine substituent with the hydrophobic pocket. Moreover, we characterized the polar intermolecular interactions of the sulfonamide group of ORY with the water molecules at the site entrance. Overall, these results suggest that ORY might be binding to the same site as PIR on the α-tubulin of T. gondii.

Models and Methods

Multiple Sequence Alignment

α-Tubulin amino acid sequences of different protozoan and nematode parasites, plants, and mammals were retrieved from the UniProt database.[44] A multiple sequence alignment was performed using T-Coffee web server.[45] BoxShade v3.21 standalone version was used for the multiple-alignment analysis and figure generation.

TgAT Model Generation

TgAT amino acid sequence (UniProt ID: P10873) was submitted to the modelling module of the SWISS-MODEL server[46] to generate a three-dimensional structure of the protein. We used the tubulin α-1B chain of Sus scrofa complexed with PIR as the template (PDB ID: 5fnv31).[47] MODELLER v9.23[48] software was used to build the C-terminal segment of TgAT (amino acids 433 to 439). The model was further subjected to a mild energy minimization with UCSF Chimera software[49] by running 100 steps of the steepest descent algorithm and 10 conjugate gradient steps. Finally, the structure assessment tool of SWISS-MODEL was employed to estimate the local and global quality of the final model.

TgAT-ORY Molecular Docking

The tridimensional structure of ORY (CID: 29393) was retrieved from the PubChem server[50] and prepared with the OpenBabel toolbox.[51] ORY was submitted to a geometry optimization and energy minimization using the MMFF94s force field[52] implemented in Avogadro package.[53] ACPYPE package[54] was used to compute ORY partial charges with the AM1-BCC method and generate the AMBER force field topology.

Rigid and Flexible Receptor Docking

ORY was docked in the PIRbs of the TgAT model using AutoDock Vina.[55] A cubic grid box of side 2.5 nm was fixed at the center of geometry of residues V4, V200, L238, S241, F255, and C316 to cover the whole binding site. The exhaustiveness of the global search was set to 50 and a maximum of 20 poses were retrieved for the docking. For the flexible receptor docking, mobility restrictions were removed to the side chains of residues L167, L238, L242, Q256, M268, C316, M318, K352, and I378. First, the flexibility of the residues was evaluated individually. Based on the resulting binding modes and scores, different combinations of the selected flexible residues were used to evaluate their contribution to the ORY interaction within the proposed binding site in TgAT. A total of 23 flexible receptor dockings were carried out in the study.

WnS Protocol

Recently, Hetényi and co-workers[56] have proposed a methodology comprising two consecutive algorithms that combines molecular docking and MD approaches to identify the most probable binding site of a ligand in a receptor. Following such methodology, first the TgAT-ORY complex, selected from the flexible docking study, was “wrapped” with a monolayer of ORY copies by performing 16 to 18 blind docking cycles with AutoDock v4.2.6.[57] The “wrapped” model was submitted to a 5 ns backbone-restrained MD simulation (MDB). We then removed the ORY copies that were displaced from their original binding site and/or presented an increase in its Lennard-Jones energy in the last simulation step. The resulting TgAT–ORY complex was further subjected to four consecutive 20 ns backbone-restrained simulated annealing MD cycles (MDBSA) to “shake” the unstable ligands and reach an elimination rate greater than 0.85. An energy minimization and NVT and NPT equilibrations were performed before each “shake” cycle. A 20 ns non-restrained conventional MD simulation for each TgAT–ORY complex was carried out, and the binding free energy of the last 5 ns was calculated with the molecular mechanics Poisson–Boltzmann surface area (MM/PBSA) method employing the g_mmpbsa(58) tool. The entropic contribution was not considered for the binding free energy calculation. All parameters used in the non-restrained simulations are described in the next section.

This protocol was performed in triplicate using the same initial TgAT–ORY complex. A clustering analysis based on the RMSD of all ORY copies in the last frame complexes was carried out to identify shared binding sites between replicates. The TgAT–ORY complexes that presented RMSD values lower than 0.5 nm between the replicates were clustered, and the MM/PBSA energy values were compared, while the complexes with only one occurrence were discarded. A total simulation time of 740 ns was yielded for the whole WnS study.

Molecular Dynamics Simulations

We carried out three independent 200 ns MD simulations of the TgAT–ORY complex employing the AMBER99SB-ILDN[59] force field implemented in the GROMACS 5.1.4 package.[60] Each replicate was simulated in a cubic box using the TIP3P water model. Sodium and chloride ions were randomly added to neutralize the charge of the system and reach a physiological concentration of 0.15 M. Amber parameters for the GTP molecule was taken from Meagher et al.[61] Each system was energy minimized employing the steepest-descent algorithm and equilibrated for 1 ns under NVT and NPT ensembles. The v-rescale coupling thermostat[62] was used to maintain the temperature at 300 K, while the pressure was isotropically fixed to 1.0 bar with the Parrinello–Rahman barostat.[63] Holonomic constraints were applied to all bonds employing the Linear Constraint Solver (LINCS) algorithm. A cut-off radius of 1.0 nm was used to compute the van der Waals and short-rage electrostatic interactions and the Particle Mesh Ewald (PME) approach was employed to approximate long-range electrostatic calculations.

For the analysis of the MD trajectories, we employed the GROMACS in-built tools to calculate the RMSD of the TgAT backbone, the RMSD of ORY heavy atoms after least-square fit to the TgAT backbone, and the TgAT backbone root-mean-square fluctuation for each simulated system. The last 190 ns of each replicate were merged into a single trajectory file to calculate the number of contacts of ORY with the TgAT residues comprising the PIRbs and the surrounding water molecules, the number of hydrogen bonds between ORY and residues N258, S241, and S237, and to carry out the RMSD clustering analysis with the GROMOS method applying a RMSD cut-off value of 0.2 nm. The MDAnalysis[64] python library was employed to calculate the distance between the center of geometry of the ORY benzene ring, the angle between the aromatic ring planes, and the number of water bridges formed at the binding site. PLIP v2.2.2[65] was used to detect the hydrophobic contacts between ORY and TgAT.

Absolute Binding Free Energy Calculations

Five representative structures of the TgAT–ORY were retrieved from the merged simulation to calculate the absolute binding free energy of the complex with the FEP method using the protocol implemented by Aldeghi et al.[66,67] In this work, electrostatic and van der Waals interactions of ORY were decoupled using a total 31 and 42 windows for the ORY-water (free) and TgAT–ORY (complex) systems. Each window was energy minimized, equilibrated for 1 ns under NVT and NPT ensembles, and subjected to 20 ns MD simulations. The free energy was computed for the last 15 ns of the ORY-water () and TgAT–ORY () simulations employing the multistate Bennett acceptance ratio (MBAR) method. ORY relative position in the PIRbs of TgAT was restrained by means of one distance (r), two angles (θ1 and θ2) and three dihedral (φ1, φ2 and φ3) harmonic potentials (kr = 1000 kcal mol–1 nm–2 and kθ,φ = 10 kcal mol–1 rad–2). The Boresch et al.[68] equation was used to analytically calculate the restraints in the decoupled ORY () system. The absolute binding free energy for each representative structure was obtained using the following equation: . A total of 0.62 and 4.2 μs were simulated for the uncoupled and coupled systems, respectively.