Homology modelling, molecular docking, and molecular dynamics simulations reveal the inhibition of Leishmania donovani dihydrofolate reductase-thymidylate synthase enzyme by Withaferin-A

Present in silico study was carried out to explore the mode of inhibition of Leishmania donovani dihydrofolate reductase-thymidylate synthase (Ld DHFR-TS) enzyme by Withaferin-A, a withanolide isolated from Withania somnifera. Withaferin-A (WA) is known for its profound multifaceted properties, but its antileishmanial activity is not well understood. The parasite’s DHFR-TS enzyme is diverse from its mammalian host and could be a potential drug target in parasites. A 3D model of Ld DHFR-TS enzyme was built and verified using Ramachandran plot and SAVES tools. The protein was docked with WA-the ligand, methotrexate (MTX)-competitive inhibitor of DHFR, and dihydrofolic acid (DHFA)-substrate for DHFR-TS. Molecular docking studies reveal that WA competes for active sites of both Hu DHFR and TS enzymes whereas it binds to a site other than active site in Ld DHFR-TS. Moreover, Lys 173 residue of DHFR-TS forms a H-bond with WA and has higher binding affinity to Ld DHFR-TS than Hu DHFR and Hu TS. The MD simulations confirmed the H-bonding interactions were stable. The binding energies of WA with Ld DHFR-TS were calculated using MM-PBSA. Homology modelling, molecular docking and MD simulations of Ld DHFR-TS revealed that WA could be a potential anti-leishmanial drug.


Introduction
Withaferin-A (WA) is among the most effective withanolide isolated from W. somnifera and has various effects like anti-bacterial, anti-inflammatory, anti-proliferative and potent anti-cancer properties [1][2][3][4]. Recently we demonstrated in vitro, that withanolides show potent anti-leishmanial activity [5] and a drastic reduction in parasite load in vivo [6].

BMC Research Notes
*Correspondence: rmusl@uohyd.ernet.in; radhemaurya@rediffmail.com 1 Department of Animal Biology, School of Life Sciences, University of Hyderabad, Prof. C.R. Rao Road, Gachibowli, Hyderabad 500046, India Full list of author information is available at the end of the article Parasite obtains folates from the host and uses its DHFR-TS and PTR1 enzymes to reduce folates to active H4 forms [22][23][24].
Hence, folate biosynthesis enzymes can be potential drug targets and molecules which inhibit any enzyme of these pathways can be a safe antileishmanial drug. Our in silico study shows that WA inhibits multiple enzymes in folate biosynthesis pathway of Leishmania parasites.

Drug-likeness prediction
Drug-likeness of WA [31,32] was calculated using molsoft server (http://molso ft.com/mprop /). A drug-likeness plot and score were obtained. Swiss target was used to predict drug target class for Withaferin A. The server, Fig. 1 Folate biosynthesis pathway, homology modelling and molecular docking: a DHFR-TS synthesizes dTMP while converting methylene THF to DHF which is converted back to THF by DHFR-TS. PTR1 converts H2 biopterin to H4 biopterin. PTR1 can reduce both pterins and folates. WA inhibits both PTR1 and DHFR-TS enzymes. b Superimposed image of the template T. cruzi DHFR-TS chain A (PDB ID: 3INV) shown in blue and modeled Ld DHFR-TS shown in green. c Substrate DHFA (red) binds to two active sites of Ld DHFR-TS where an electrostatic channel is formed and substrate channeling between both the active sites is observed. Competitive inhibitor MTX (blue) competes with DHFA (red) and binds to two active sites of Ld DHFR-TS. Inhibitor WA (yellow) is binding to Ld DHFR-TS enzyme by blocking the electrostatic channel. d Substrate DHFA (red), Competitive inhibitor MTX (blue) and inhibitor WA (yellow) binding to Hu DHFR enzyme. e Substrate dUMP (red), inhibitors MTX (blue) and WA (yellow) binding to Hu TS enzyme using a combination of 2D and 3D similarity measures, compares the query molecule to a library of 280,000 compounds active on more than 2000 targets of five different organisms [33].

Molecular dynamic simulation
MD simulation of Ld DHFR-TS, Hu DHFR and Hu TS, and their WA complexes were performed in Gromacs 5.0. (http://www.groma cs.org/) [34]. The topological parameter of the ligand was obtained from ATB server (https ://atb.uq.edu.au/) [35]. Initially, protein or its complex was kept in a cubic box filled with water using SPC/E water models. The system was energy minimized using GROMOS54a7 force field [36] and equilibrated at 300 K using V-rescale for 200 ps as NVT ensemble followed by equilibration at 1 atm pressure using Parrinello-Rahman algorithm as NPT ensemble for 200 ps. The equilibrated conformation was further extended for production simulation for 25 ns. LINCS algorithm was applied for bond constraints with distance cut-off using Verlet during simulation. Root mean square deviations of atomic coordinates during the simulation from their respective initial coordinates were calculated using the gmx_rms tool in Gromacs and binding energies were calculated using MM-PBSA [37].

Sequence alignment and homology modeling
The sequence similarity between Hu DHFR and Ld DHFR-TS was found to be 25.13%, and between Hu TS and Ld DHFR-TS, it was 54.63% suggesting that Ld DHFR-TS could be a valid drug target (Additional file 1: Figs. S2, S3). The amino acid sequence of Ld DHFR-TS was blasted against PDB-BLAST database for identifying an appropriate template for homology modeling. T.cruzi DHFR-TS showed 67.32% identity with the target protein and was selected as a template (Additional file 1: Fig.  S4). Quality of the model generated by Swiss-model was verified using different tools (Fig. 1b) (Additional file 1: Table S1). The selected model showed 0.2% of residues in disallowed regions of Ramachandran plot (Additional file 1: Fig. S5, Table S2) with GMQE score of 0.82 and QMEAN score of − 2.25 (Additional file 1: Fig. S6).
The generated model is a homo-dimer protein of α + β class. The protein consists of 4β-sheets, 3βαβ units, 5β-hairpins, 19β-strands, 21α-helices (Additional file 1: Fig. S7). Similar numbers of secondary structural elements were found in T. cruzi DHFR-TS and RMSD between the template and generated model was calculated to be 0.625 Å.

Molecular docking studies
To know the active site of Ld DHFR-TS, it was first docked with its substrate DHFA and found that it has two active sites, one in DHFR and other in TS domain. TS active site is located 40 Å away from DHFR active site [38][39][40].  (Fig. 1c).
Crystal structure of Hu DHFR (4m6k) was docked with WA and was superimposed with Hu DHFR ternary crystal complex of MTX and NADPH (1u72) and crystal structure of Hu DHFR complex of NADP+ and folate (4m6k). The results showed all three ligands viz. WA, DHFA, and MTX are binding in the same pocket. The ligand WA also competes for the active site and might be acting as a competitive inhibitor. The binding energy of WA is − 41.4 kJ/mol (Fig. 1d).
Crystal structure of Hu TS (1hzw) was docked with WA and later superimposed with Hu TS complex of dUMP and MTX (5x66). We observed that WA is binding at the same site like MTX. The residues Phe 80, His 196, Leu 221 and Asn 226 were forming H-bonds with WA and binding energy of WA was − 39.8 kJ/mol. The ligand WA was again competing for the active site and might be acting as a competitive inhibitor (Fig. 1e) The binding energy of Ld DHFR-TS with WA is higher than Hu DHFR and TS. Moreover, WA could be a better drug than MTX because of its high binding energy.

Molecular dynamic simulations of enzyme-inhibitor complexes
To characterize the stabilizing interactions and to evaluate binding energies of WA with Ld DHFR-TS, Hu DHFR and Hu TS, MD simulation of proteins and protein-WA complexes were carried out. The analysis of root mean square deviations (RMSD) showed all proteins attained almost stable conformations (Fig. 2a-c) with comparable RMSD values. Addition of WA did not show much change in RMSD of Hu DHFR whereas RMSD of Ld DHFR-TS slightly increased. RMSD of ligand alone was around 0.15 nm in all proteins suggesting that bound conformation was stable. Further, root mean square fluctuations (RMSF) of individual residues were calculated by considering their Cα atoms as a reference (Fig. 2df ). RMSF of β5-loop in DHFR domain and β1′ and β4′ For further quantitative binding, energies of ligand were calculated by MM-PBSA using the last 10 ns of simulation data where RMSD of proteins were found to be more stable ( Table 1). The analysis indicates that binding affinity of WA is more towards Ld DHFR-TS than Hu DHFR or Hu TS.

Discussion
Interestingly, Leishmania dhfrts− mutants are unable to survive in mammalian host [41]. Deletion of PTR1 gene is lethal in promastigotes, indicating an essential role for unconjugated pteridines [20][21][22][23]42]. PTR1 expression provides a potential 'metabolic by-pass' of DHFR-TS inhibition and allows a partial or complete reversal of anti-pteridine inhibition in the promastigote stage of parasites [20,21]. PTR1 activity in L. major promastigotes is lower than in L. donovani and L. mexicana. L. major is more sensitive to MTX suggesting the role of PTR1 as a metabolic-bypass in L. donovani and L. mexicana [18,19]. 3D structures of DHFR-TS and PTR1 of parasite and Hu DHFR have provided a strong base to design new inhibitors which are selective for parasite alone [43,44].
Recently, we reported that WA inhibits Ld PTR1 enzyme activity and molecular docking studies of WA showed high binding affinity with PTR1. Enzyme assay with purified PTR-1 revealed that WA inhibits enzyme activity through uncompetitive mode [45]. The present molecular docking study reveals that the binding energy of WA with Ld DHFR-TS is higher than Hu DHFR, Hu TS enzymes and WA inhibits Ld DHFR-TS same as the PTR-1 enzyme. Thus it could be concluded that binding affinity of WA with multiple enzymes (DHFR-TS and PTR1) of folate biosynthesis pathway of parasites could make WA an effective anti-leishmanial drug.

Limitation
Due to the lack of purified DHFR-TS enzyme, the current study could not include enzyme assay. However, enzyme assayed from parasite lysate with WA has shown the inhibition activity reported earlier [45].  Figure S6. Local quality estimate of (A) modelled Ld DHFR-TS and (B) reference T.cruzi DHFR-TS obtained from Swiss model. Figure S7. Secondary structures of (A) modelled Ld DHFR-TS and (B) reference T.cruzi DHFR-TS obtained from PDB sum. Table S1. Features of the generated Ld DHFR-TS model from Swiss model. Table S2. Ramachandran plot Statistics from PROCHECK results for modelled Ld DHFR-TS protein and reference T. cruzi DHFR-TS protein. Table S3. Drug likeness properties of WA from molsoft.