Febrifugine analogues as Leishmania donovani trypanothione reductase inhibitors: binding energy analysis assisted by molecular docking, ADMET and molecular dynamics simulation

Visceral leishmaniasis affects people from 70 countries worldwide, mostly from Indian, African and south American continent. The increasing resistance to antimonial, miltefosine and frequent toxicity of amphotericin B drives an urgent need to develop an antileishmanial drug with excellent efficacy and safety profile. In this study we have docked series of febrifugine analogues (n = 8813) against trypanothione reductase in three sequential docking modes. Extra precision docking resulted into 108 ligands showing better docking score as compared to two reference ligand. Furthermore, 108 febrifugine analogues and reference inhibitor clomipramine were subjected to ADMET, QikProp and molecular mechanics, the generalized born model and solvent accessibility study to ensure the toxicity caused by compounds and binding-free energy, respectively. Two best ligands (FFG7 and FFG2) qualifying above screening parameters were further subjected to molecular dynamics simulation. Conducting these studies, here we confirmed that 6-chloro-3-[3-(3-hydroxy-2-piperidyl)-2-oxo-propyl]-7-(4-pyridyl) quinazolin-4-one can be potential drug candidate to fight against Leishmania donovani parasites.


Introduction
Leishmaniases are geographically widespread, neglected tropical diseases which affect about 350 million people from total 98 countries and three territories on five continents (Alvar et al., 2012). The disease is caused by more than 20 species of leishmania and is transmitted by the bite of infected female sandfly vector of Phlebotomine genera (Chappuis et al., 2007). The forms of leishmania which spend its life in female sandfly are extracellular flagellated promastigotes, while it replicates as intracellular aflagellated amastigotes in mononuclear phagocytes of the mammalian host (Killick-Kendrick, 1990). Depending on the leishmania species, the disease symptoms may range from self-healing cutaneous lesions to the metastating muco-cutaneous form and to the more fatal visceral form known as kala-azar or visceral leishmaniasis (VL). Kala-azar is the most severe form of the leishmaniasis and accounts for 58,200 new cases and over 30,000 death per year, while 220,000 cases were reported for cutaneous leishmaniasis annually (Alvar et al., 2012). In the Indian subcontinent and East Africa, the primary cause of VL infection is Leishmania donovani while in Mediterranean basin, Middle East and New World, VL infection is caused by Leishmania infantum. The disease symptoms include fever, splenomegaly, hepatomegaly, cachexia, epistaxis and pancytopenia. Majority of VL cases are found in India, Bangladesh, Brazil, Sudan, South Sudan and Ethiopia (Alvar et al., 2012). The situation of VL is especially severe in Bihar state in eastern India where some districts have faced the worst epidemic since the 1970s (Prajapati, Mehrotra, Gautam, Rai, & Sundar, 2012).
Antileishmanial drugs, currently in use are either expensive or require long-term treatment and most of them are associated with adverse events. Pentavalent antimonials have been used as primary drug for the treatment of leishmaniasis worldwide. However, in India, reports from Bihar suggested high level of resistance with only one-third patients responding to the treatment. Thus, antimonials are not in longer use in the endemic zone of India (Freitas-Junior, Chatelain, Kim, & Siqueira-Neto, 2012;Perry et al., 2011). The efficacy of amphotericin B and liposomal amphotericin B is high (>95%), but amongst these two drugs, amphotericin B has shown significant nephrotoxicity while its liposomal formulation causes rigours and chill during infusion (Freitas-Junior et al., 2012). Paromomycin, an aminoglycoside antibiotic isolated from Streptomyces rimosus, possesses antiprotozoal activity and is even effective against leishmania infection. It had shown cure rate of 94.6% when injected intramuscularly in patient suffering from VL (Sundar, Jha, Thakur, Sinha, & Bhattacharya, 2007). The major advantage of this drug is its cost-effectiveness, but requirement for intramuscular injections for three weeks, monitoring of serum transaminase and lack of adequate data regarding its safety in pregnancy are its drawbacks. There is also an increased resistance to newly introduced antileishmanial oral drug namely miltefosine. It has shown 7-10% patients relapse at 6 month and up to 20% at 12 month in the clinical trials conducted in the Indian subcontinent Rijal et al., 2013). In the recent studies, clinical isolates of L. donovani have shown high fitness and tolerance against miltefosine drug which also supports resistance phenomenon (Prajapati et al., 2013). Toxicity, high cost and increasing resistance of antileishmanial drugs entails development of new drug against leishmania with good efficacy and safety profile (Sundar & Chakravarty, 2014).
Trypanothione reductase (TR) (PDB ID-2JK6), a central enzyme involve in the redox metabolism of leishmania parasite, is a 511 amino acid residues containing NADPH dependent flavoprotein (Fairlamb & Cerami, 1992). Flavin adenine dinucleotide (FAD) is a cofactor which tightly bound to each subunit of TR dimer and help in transfer of electron during catalysis (Baiocco, Colotti, Franceschini, & Ilari, 2009;Fairlamb & Cerami, 1992). There are two active sites occupied by TR dimer. The first active site is NADPH binding site. Apart from this, there is a presence of second active site at the interface of chain A and chain B. TR helps in maintenance of trypanothione in its reduced state (TSH 2 ) by the transfer of electron from NADPH to oxidized trypanothione through FAD and cysteine residue (Cys-52 and Cys-57) (Fairlamb, Blackburn, Ulrich, Chait, & Cerami, 1985;Meiering et al., 2005). Cys-52 and Cys-57 are highly conserved amino acid residue present in the active site of TR and they play important role during catalytic action of TR (S. Rai, Dwivedi, & Goyal, 2009). The absence of TR in human and its important role in parasite antioxidant system makes it a novel drug target for drug discovery against leishmania parasites (Meiering et al., 2005).
Virtual screening of a series of compounds has become an initial approach to reduce the virtual space of chemical ligand up to a manageable level where we could synthesize and check the biological effectiveness of selected compound (Seeliger & de Groot, 2010). Structure-based virtual screening assist the protein interactions with ligand at the atomic level and it allow us to identify the behaviour of chemical molecules in the binding pocket of target protein (Hazra et al., 2013;Verma et al., 2012). Molecular mechanics, the generalized born model and solvent accessibility (MM-GBSA) is an additional module of Schrödinger suite, used as a post-docking scoring protocol to correctly rank the potent inhibitor molecule against target protein (Lyne, Lamb, & Saeh, 2006). MM-GBSA uses molecular mechanics, the generalized Born model and Solvent accessibility method to obtain free energy from structural information to determine the relative binding free energies (ΔG bind ) in bimolecular complexes. ADMET predictor (ADMET predictor TM V7.2) is used for the advanced predictive modelling of ADMET properties and simultaneously to simulate the absorption profile of drug (Okumu, DiMaso, & Löbenberg, 2009). The mathematical model in the programme uses molecular descriptor value to generate estimate for each ADMET properties (Fraczkiewicz et al., 2015;Singh, Gupta, Singh, Singh, & Misra, 2013). Schrödinger QikProp is used to predict pharmaceutically applicable properties (ADME) either of single compound or compounds in batches (Ioakimidis, Thoukydidis, Mirza, Naeem, & Reynisson, 2008). AutoDock4.2 is another automated and robust docking algorithm based on the Lamarckian Genetic Algorithm (LGA) with significant success rates in the virtual screening of ligand. It also has free energy scoring function which uses AMBER force field to estimate the binding energy of ligand to its target (Morris et al., 2009;Rashid, Kuppa, Kunwar, & Panda, 2015;Venghateri, Gupta, Verma, Kunwar, & Panda, 2013).
Febrifugine is an alkaloid drug mainly obtained from Dichroa febrifuga (Chinese herb), a saxifragaceous plant. It is used for the treatment of malaria fever (causal organism: Plasmodium) for more than 2000 years in Chinese traditional medicine, and no parasitic resistance has been reported against the febrifugine obtained from D. febrifuga (Takaya et al., 1999). By conducting in vivo studies, researchers have demonstrated the less toxic and better therapeutic behaviour of febrifugine analogues, which supports the potential animalarial activity of febrifugine analogues (Zhu et al., 2010). As febrifugine analogues were never tested for the treatment of VL infection, we have tried to verify its efficacy against TR to kill L. donovani protozoan parasite and control VL infection. In this article, we have emphasized the interaction of febrifugine analogues and TR enzyme of L. donovani antioxidant system. Structure-based virtual screening was further followed by toxicity determination, ADME properties prediction and binding free energy calculation using ADMET predictor, QikProp and MM-GBSA module of Schrödinger suite, respectively. The robustness of docking protocol was analysed through ROC curve and BEDROC analysis using enrichment calculator of Schrödinger suite. This study also includes docking validation for best two analogues and reference inhibitors using another docking module "AutoDock 4.2". Furthermore, we have employed molec-ular dynamics (MD) using Desmond to expand better understanding regarding interaction between proteins and ligand ( Figure 1).

Molecular docking
To get potential antileishmanial drug, docking scores were calculated using Glide (Glide) module of Maestro 10.2 (Maestro). In this study, initially the compounds with greater than and equal to 80% similarities to febrifugine were extracted from PubChem database in SDF format and were subsequently used as the library of febrifugine analogues. Then the docking was performed for these analogues against TR of Leishmania parasites with a variety of sequential computational phases.

Protein preparation
The typical protein structure imported from PBD database is not suitable for binding affinity calculation because they may consist of co-crystallized ligands, heavy water molecule, metal ions and co factors. PDB structures also may have missing atoms, missing loops and connectivity information which must been assigned along with bond order and formal charges. In present study, the structure of trypanothione reductase (PDB ID-2JK6) was obtained from PDB and prepared using protein preparation wizard (PrepWizard) in Maestro (Protein Preparation Wizard 2015-2; Epik version 2.4, Impact version 5.9, & LLC). Initially H-bond was optimized and then the whole protein structure was allowed to relax using Impref module of Impact (Impact) and finally the receptor protein was optimized by applying OPLS_2005 force field.

Grid generation
The grids for molecular docking were made using Glide (Glide) to decide a site on the target protein where ligand has to interact during molecular docking process. There are two active sites present in TR, one active site is occupied by FAD and another active site is at the interface of chain A and B of TR dimer. In present study, grid was made at the interface of AB chains dimer. The amino acid residues of active site involves in bonding or non-bonding interactions in the respective chain are Val-58, Lys-61, Leu-62 and Thr-65 from the B chain and Phe-396, Thr-397, Pro-398, Leu-399, Met-400, His-461, Pro-462, Thr-463 and Ser-464 from A chain (Verma et al., 2012). All aforementioned amino acid residues surround active sites Cys-52 and Cys-57 residue of B chain to form the binding pocket (Rai et al., 2009). The size of grid was kept default at 20 Å.

Docking analysis
Docking was conducted using Glide v6.7 (Glide) in three steps viz. High throughput virtual screening (HTVS), Standard precision (SP) and extra precision (XP) mode to screen the compound with highest binding energy against target protein. Amongst these three docking mode HTVS and SP uses self-same scoring function where former minimizes the number of intermediate conformation throughout the docking and also reduces the thoroughness of the endmost torsional refinement and sampling. XP also uses more sophisticated scoring function than SP to perform more extensive sampling than SP. Total 23,492 stereoisomers of febrifugine analogues including two reference ligands were docked against TR by HTVS mode of Glide to superficially select the best interacting ligand. Later on, top 471 ligands were selected and put forward for SP Glide docking analysis. The top 109 ligand showed best negative docking score for SP were introduced for XP docking along with two reference ligands (previously mentioned).

ADME properties
Schrӧdinger QikProp v4.4 (QikProp) was used to calculate ADME properties for 108 febrifugine ligands and clomipramine. QikProp module uses ligand either from project table or uses the LigPrep files saved on computer storage while the output result consist of a number of principal descriptor and ADME prediction.

Molecular mechanics, the generalized born model and solvent accessibility
MM-GBSA (Prime) another module of Schrödinger software was used to obtain relative binding energy of selected ligand from structural information in bio molecular complexes. During this study, the active site of protein was set to adjust itself up to a distance of 5 Å for ligand accordingly. Total 108 ligands and clomipramine, were subjected to the MM-GBSA analysis. This operation import XP docking "out.maegz" file of protein and ligand complex, which results in ranking of ligand based on calculated binding energies (MM-GBSA DG binds). Amongst aforementioned ligands, top 29 ligands showing good binding energy were subjected to ROC curve plot and enrichment study. Furthermore, amongst 29 ligands, the best two ligands qualifying all aforesaid screening parameters (docking, toxicity and binding energy) were subjected to another docking protocol viz. AutoDock4.2 MD simulation using Desmond.

Enrichment calculation and validation of docking protocol
The robustness and ability of docking protocol to retrieve active compounds amongst a set of inactives (Decoys) is based on enrichment calculation and receiver operator characteristics (ROC) analysis. The top 29 febrifugine analogues, showing best docking score with less toxicity were selected as active ligand. A total of 36 decoys per active ligands were prepared using DecoyFinder (Cereto-Massague et al., 2012). Later on, all 29 active febrifugine analogues were mixed with 1044 decoy compounds (Zinc library) and docked against TR dimer using SP and XP mode of Glide. Separate ROC curve were plotted for both SP and XP docking, using enrichment calculator of Schrӧdinger Inc. The value of enrichment factor and BEDROC is also incorporated in this study. We used α = 160.9, α = 20.0 and α = 8.0 for the BEDROC calculation.

Molecular docking validation using AutoDock 4.2
To compare the performance of docking technique, another docking programme namely AutoDock 4.2 was used. The three dimensional coordinates of the ligand FFG2 (NCI ID: 21051082) and FFG7 (NCI ID: 76799884) along with Clomipramine and 4-(4,4,8-Trimethyl-7-oxo-3-oxabicyclo[3.3.1] non-2-yl)-benzoic acid methyl ester were generated from the 2D SDF file using PyMol software, followed by the energy minimization and optimization using the "Prepare Ligand" in the AutoDock4.2 for the docking studies. The X-ray crystallographic structure of TR protein from Leishmania infantum (PDB ID: 2JK6) resolved at 2.95 Å was retrieved from the Protein Data Bank (Baiocco et al., 2009). In order to identify putative binding site of febrifugine inhibitor on TR dimer, water molecules and SO 4− ions were removed from the crystal structure. The hydrogen atoms and Gasteiger charges were added to TR protein using UCSF Chimera software (Pettersen et al., 2004) and Autodock4.2, respectively (Morris et al., 2009).
Here, we employed the blind and local docking protocol, to identify the putative binding site of ligand, as used in previous docking studies of anticancer drugs with tubulin dimer Rashid et al., 2015;Venghateri et al., 2013). For blind docking (Hetényi & van der Spoel, 2006;Iorga, Herlem, Barré, & Guillou, 2006) the TR protein dimer was enclosed in a grid box of size 126 Å × 126 Å × 126 Å with a spacing of 0.7000 Å. The entire TR dimer was taken as rigid, while febrifugine ligands were taken as flexible molecules. Further the centre position of grid box was set at the interface TR protein.
During molecular docking simulation, the LGA was used with the default parameters; g_eval set to 2,500,000 (Medium). Total of 5 independent flexible ligand docking jobs were conducted, each with 100 LGA runs, which produced 500 output conformations in each case. The febrifugine ligand along with reference ligands were placed at random locations for each docking simulation run. The febrifugine and reference ligands bind at interface of dimer of chain A and chain B as similar to chromene-2-thione analogues and other TR inhibitor (Venkatesan & Dubey, 2012;Verma et al., 2012) For local docking, the grid box was set to cover the active site along with maximum search space around the shared binding pocket. Here, a grid box with a size of 126 Å × 126 Å × 126 Å with a spacing of .375 Å was constructed around the active site region. The LGA was also employed with the default parameters; g_eval set to 2,500,000 (Medium). Total 50 independent flexible ligand docking jobs were conducted each of 100 LGA runs for both FFG2 and FFG7 ligand along with reference ligand, which produced a total of 5000 output conformations for each.
The Autodock4.2 (Morris et al., 2009) was used to perform clustering of binding conformations of ligand based on the scoring functions with information about the frequency of occurrence, mean energies, inhibition constant and root mean square deviations within the cluster. The resulted 5000 output conformations were further clustered using an all-atom RMSD cut-off of 4 Å. The clusters having more than 100 conformations (cutoff number) were selected and analysed further. The clusters were then compared on the basis of the cluster size and binding energy calculated by AutoDock 4.2 scoring function. The TR protein and febrifugine ligand structural images and hydrogen bonding interactions analysis were performed using PyMol (DeLano, 2002) and UCSF Chimera (Pettersen et al., 2004).

MD simulation
MD simulation is widely used method to understand the properties of the structures and the microscopic interactions between them (Alder & Wainwright, 1959). Fullatom MD simulations were performed for each system using Desmond (Desmond Molecular Dynamics System, 2015. Maestro-Desmond Interoperability Tools, 2015). The whole process completed in three stages viz. system builder, minimization and MD. Solvation builder was used to select SPC as a solvent model with orthorhombic boundary box and charges were neutralized using sodium and chloride ions. Before simulation, we ensured that the system has no steric clashes by subjecting it to process of energy minimization. The system equilibration was done via NVT and NPT ensemble using shake algorithm and by bringing the temperature up to 300 K and pressure up to 1 bar. After equilibration 10 ns MD simulation was executed for the protein and protein-ligand complex. The bond length was constrained using LINCS algorithm (Hess, Bekker, Berendsen, & Fraaije, 1997). The simulation was carried for 10 ns and trajectories were analysed. After completion of MD simulation, the resulting md.out.cms file was used to calculate RMSD of protein-ligand complex, Protein RMSF, Ligand RMSF, protein ligand interaction and ligand torsion study.

Toxicity analysis using ADMET predictor
Amongst selected toxicity parameters, MRTD represent a qualitative assessment and shows the maximum recommended therapeutic dose administered as an oral dose. Any effective oral drug must have the MRTD value greater than 3.16, which confirms that the selected compound is "inactive" with fewer side effects and vice versa (Matthews, Kruhlak, Benz, & Contrera, 2004;Singh et al., 2013). Amongst the finalized 29 ligand, total 13 compounds including clomipramine showed MRTD value greater than 3.16 while the remaining has the reverse value. The second parameter viz. hERG (The human Ether-à-go-go-Related Gene) represent the gene which encodes potassium ion channel, causes normal repolarization of the cardiac action potential. The blockage of this channel leads to cardiac arrhythmias (Sanguinetti & Tristani-Firouzi, 2006). The IC 50 value of less than or equal to 10 μmol/L were toxic, while reverse is non-toxic and unable to block hERG channel (Sanguinetti & Tristani-Firouzi, 2006). Amongst febrifugine analogues, FFG3 has shown the highest hERG value followed by FFG1, FFG2, FFG20, FFG7 and FFG17. Whereas, clomipramine has shown the hERG value greater than FFG3. CARB shows the abbreviated representation of chromosomal aberration in response to chemical compound. When febrifugine analogues were evaluated for this parameter FFG1, FFG2, FFG4, FFG13, FFG14, FFG20 and FFG24 were found to be non-toxic against clomipramine. Phospholipidosis (PHOS) is a lysosomal storage disorder, causes accumulation of phospholipid which results in impaired neuronal signalling. Total 25 ligands were found to be non-toxic, except clomipramine.

Principal descriptor and ADME analysis using QikProp
Total 17 ADME properties for febrifugine analogues were calculated, using Schrӧdinger QikProp module to validate safety profile of compounds (Table 3) (Butina, Segall, & Frankcombe, 2002;Yamashita & Hashida, 2004). The ADME properties include, predicted polarizability in cubic angstrom (PolrZ), predicted hexadecane/gas partition coefficient (Log PC16), predicted octanol/gas partition coefficient (Log Poct), predicted water/gas partition coefficient (Log Pw), predicted octanol/water partition coefficient (Log Po/W), predicted aqueous solubility (Log S), conformation independent predicted aqueous solubility (CI logS). Other included ADME properties were Caco-2 permeability in nm/s (PCaco), predicted brain/blood partition coefficient (Log BB), predicted skin permeability (Log Kp), number of likely metabolic reaction (Metab) and prediction of binding to human serum albumin (Log Khsa). Amongst ADME properties, predicted qualitative human oral absorption (HOA, value ranges from 1 for low, 2 for medium and 3 for high), predicted human absorption on percentage scale (PHOA), number of violation of  . Whereas, total solvent-accessible volume in cubic angstroms using a probe with a 1.4 Å radius (Mol vol), estimated number of hydrogen bonds that would be donated by the solute to water molecule in an aqueous solution (Donor HB), estimated number of hydrogen bonds that would be accepted by the solute from water molecule in an aqueous solution (accpt HB) and globularity descriptor (glob) and van der Waals surface area of polar nitrogen and oxygen atom and carbonyl carbon atoms were also fulfil the principle descriptor properties. All compounds have shown positive result for all principle descriptor along with clomipramine except FFG20 (Boelaert et al., 2009;Ioakimidis et al., 2008;Singh et al., 2013) (Table 4).

MM-GBSA of febrifugine analogues
The binding energy of top 29 febrifugine analogues and clomipramine were calculated from their structural information in binding pocket of biomolecular complexes using Prime MM-GBSA and it was noted that FFG7and FFG2 showed highest negative MM-GBSA dg bind value of −44.5 and −40.3, respectively. On the other side, clomipramine showed higher MM-GBSA dg bind value (−21) as compared to aforesaid febrifugine analogues (Table 1).

ROC curve plot and enrichment study
The ROC curve was plotted using the pose viewer (pv.maegz) file of docking output (29 actives and 1044 decoys) of SP and XP result (Figure 2). The ROC curve plot analysis resulted into better performance for XP  −2.5/6.5 −6.5/.5 −6.5/.5 25 to >500   Figure 4(A)-(B). The binding energies of FFG2 and FFG7 were estimated to be -12.16 kcal/mol and -11.42 kcal/mol, respectively. While the binding energies for clomipramine and 4-(4,4,8-Trimethyl-7-oxo-3-oxabicyclo[3.3.1]non-2-yl)-benzoic acid methyl ester were found to be −11.29 and −8.64, respectively by docking analysis. The binding pocket of TR protein around 4 Å distance of FFG2 is amphipathic in nature, i.e. contains both hydrophobic as well as hydrophilic residues (Figure 3(A); Supplementary Table SI1). These residues are involved in bonding and nonbonding interactions with FFG2 are as follows: Lys-61, Pro-336, Ile-339, Asn-340, Glu-436 from B chain and   Table S1). The FFG2 is stabilized by hydrogen bonding with Lys-61(2.79 Å) and Pro-336 (2.16 Å) from B chain and Thr-463(2.44 Å), Glu-466(2.10 Å) and Glu-467 (2.16 Å) from A chain (Figure 3(A), Supplementary  Table SI2). It has observed that the FFG2 is surrounded and interact mostly with the chain A residues and few residues with chain B at the interface. The FFG2 shows hydrogen bonding with residues Glu-466 and Glu-467, these residues form the hydrogen bonds with the active site His-461 of B chain which forms the core of the hybrid transfer region. The interactions of FFG2 with these residues (Figure 3  might be changing the active site structure, which ultimately leads to the inactivation of the TR enzyme. Further, we also carried out the docking of FFG7 inhibitor with TR enzyme using similar approach. Docking analysis suggested that FFG7 prefer similar binding site at the interface of dimer on TR protein (Figure 4(A) and (B)) as for FFG2 (Figure 4(A) and (B)). The estimated binding energy of FFG7 was found to be −11.42kcal/mol. The binding pocket of FFG7 also consists of both hydrophobic and hydrophilic amino acids around the 4 Å distance which are as follows:   Table S1 inline). The FFG7 is stabilized by hydrogen bonding interactions with as Glu-466 (1.63 Å), Glu-467 (2.10 Å) and Ser-470 (1.94 Å) from A chain (Figure 3(B), Supplementary Table SI2), it doesn't show any interactions with the residues of chain B.

from A chain (Supplementary
Further, the electrostatic potential was calculated over the docked complex of TR dimer and FFG2 and FFG7 (Figure 4(C) and (D)), to understand the binding mode of ligand at the interface of TR protein using PyMol. The red, blue and white colour shows the electro-negative, electro-positive and neutral potential, respectively (Figure 4(C) and (D)). The electrostatic potential surface of TR dimer shows that, the FFG2 and FFG7 bound at the interface cavity formed between chain A and chain B (Figure 4(C) and (D)) as similar to antileishmanial chromene-2-thione analogues (Verma et al., 2012). Here, both the FFG2 and FFG7 are docked at the same binding pocket and show the similar interactions with the residues at the binding site.

MD simulation study
MD simulation study was conducted to verify the interaction between TR and best two febrifugine analogues qualifying all screening parameters using Desmond v2014.2. In the beginning phase both the complexes viz. FFG2 and FFG7 were subjected to system builder separately. SPC was selected as a water model and sodium or chloride ions were added for charge minimization. FFG2 model system contained total 102,089 atoms and 29,008 water molecules. While the FFG7 model system consist of 102,212 atoms and 29,047 water molecules. The overall FFG2 complex was neutralized using 109 NA + and 89 Cl − ions while that of FFG7 using 105 Na + and 81 Cl − ions. Further steps were followed by Broyden-Fletcher-Goldfarb-Shanno algorithm. MD simulations for both these complexes were performed at 300°K and 1.0 bar. The structural stability of trypanothione reductase (PDB ID: 2JK6) and the protein ligand complex (FFG2 and FFG7) were examined during simulation by analysing root mean square deviation (RMSD). RMSD of protein during the simulation measure the stability of system as it gives an idea about type of binding of the ligand and the structural changes that occurred during simulation. The RMSD of FFG2 and FFG7 protein ligand complex for 10 ns is around 1.6 and 2.0 Å, respectively, which is within the satisfactory range of 0-3 Å. But if we consider comparative stability then FFG2 complex has shown more stability characteristic than FFG7 complex (Figure 5(A) and (B)).
Root mean square fluctuation (RMSF) is another parameter which gives an idea about the conformational and structural changes occurs along the protein side chain. The RMSF of protein has shown flexibility within a range of .5-2 Å, shows that the binding of FFG2 leads to stable affinity for protein active site. While the amino acid residue in protein side chain (480th amino acid residue) has shown fluctuation up to 4 Å due to its terminal position. The N and C terminal residue have most of the fluctuation while the central core region is stabilized. The value of RMSF was also calculated for each residue. Residues Gly481, Val 441and Gly442 have shown maximum fluctuation during simulation as they lie near active site of protein. The RMSF of protein complexed with FFG7 has shown different results. The protein has shown RMSF value within a flexible range from 0.6 to 4.0 Å and it shows higher flexibility than FFG2. This result shows that binding of FFG7 to protein active site leads to moderate stability. The amino acid residue have shown maximum fluctuation are Gly300, Val301, Val460 (Figure 5(C) and (D)).
Ligand RMSF is another parameter to characterize changes in ligand position. The protein-ligand complex was first aligned on the protein backbone and then the ligand RMSA was measured on the ligand heavy atom. In FFG7 only second atom has shown fluctuation up to 2 Å while remaining atom has moderate fluctuation. While Atom no 11th, 21st, 25th and 31st of FFG2 ligand has shown a major fluctuation up to 2 Å. (Supplementary Figure SI 1A-B).
The protein ligand interaction was also monitored throughout the simulation study and summarized for FFG2 and FFG7. Protein ligand interaction were categorized into hydrogen bond, hydrophobic bond, ionic and water bridges. Amongst all these interaction H-bonding is important because of their strong influence on drug specificity, metabolization and adsorption. For FFG7, Glu466 made hydrogen bond, ionic and water bridges (Supplementary Figure SI2) beside this Ser470 and Arg472 have also shown hydrogen bond and ionic interaction. Amino acid residue Glu466 and Glu467 has shown to be participated in the formation of hydrogen bond and water bridges while Glu466 has also shown ionic interaction with FFG2 (Supplementary Figure SI3).

Discussion
In the Indian subcontinent, the treatment option strongly depends upon pentavalent antimonials, amphotericin B deoxycolate (dAmB), paramomycin and miltefosine . All these drugs are either expensive or require long-term treatment and also associated with resistance to leishmanial parasite (Chakravarty & Sundar, 2010;Perry et al., 2011;Prajapati et al., 2011Prajapati et al., , 2012. Liposomal amphotericin B has shown good efficacy (more than 95% in a single dose) and less toxicity, but its high cost and requirement of cold chain, pose serious problem for its development in the field (Mukhopadhyay et al., 2011). Febrifugine is an active antimalarial drug in Chinese tradition for more than 2000 years, which was isolated from Chinese herb chang shan (D. febrifuga) (Jiang et al., 2005). Febrifugine derivatives have shown wide range of activity to treat cancer, fibrosis, malaria and inflammatory diseases. Amongst the febrifugine derivative, halofugine has shown inhibition effect against development of Th17-driven autoimmunity in mouse model of multiple sclerosis (Keller et al., 2012). Many natural compounds and their chemically synthesized derivatives have shown antileishmanial activity by targeting TR enzyme present in L. donovani parasites. Plumbagin is one of them, has shown strong inhibition activity against leishmania parasite (Sharma, Shukla, Das, & Dubey, 2012). Betulin is a cost-effective natural product, having antileishmanial properties even at very less concentration. It induces apoptosis by stimulating reactive oxygen species generation against leishmania promastigotes leads to cure (Saudagar & Dubey, 2014). Therefore, above argument enforce the urgent need of antileishmanial drugs having main skeleton of natural compounds to avoid associated side effects during treatment. By keeping these abovementioned parameters, we have screened febrifugine analogues against trypanothione reductase to verify their inhibitory activity to kill L. donovani parasites.
The present study focused on the virtual screening, toxicity analysis and MD simulation analysis to find novel L. donovani parasites inhibitor. Virtual screening was used to predict the interaction of ligand within the binding pocket of target enzyme which results in docking score. Higher the negative value of docking score, stabilizes the association and proportionally the parasite clearance (Sahoo et al., 2014;Verma et al., 2012). In this study, we evaluated the relative binding affinity of 8813 febrifugine analogues against binding pocket of TR enzyme of L. donovani parasites. Initially, the evaluation of binding affinity was carried out by molecular docking approach in two sequential docking mode (HTVS and SP) and finally result was concluded by XP mode of docking for 109 ligands. We also evaluate the same for clomipramine (wellknown TR inhibitor) and 4-(4,4,8-Trimethyl-7-oxo-3oxabicyclo[3.3.1]non-2-yl)-benzoic acid methyl ester for experimental validation. When the results were analysed, it concluded that total 108 febrifugine analogues have higher docking score than clomipramine and 4-(4,4,8-Trimethyl-7-oxo-3-oxabicyclo[3.3.1]non-2yl)-benzoic acid methyl ester. Amongst 108 febrifugine analogues the best interacting ligand was 6-[[3,4,5-trihydroxy-6-(hydroxymethyl)tetrahydropyran-2-yl]amino]-1H-quinazolin-4-one having docking score -7.96. But when all these best interacting compounds (total 108 febrifugine analogues and clomipramine) were evaluated on toxicity parameter using ADMET predictor and QikProp, we noticed that only 10 compounds were qualifying all toxicity parameter, while remaining have shown toxicity on few parameters and many analogues have shown high toxicity. Later on, we finalized total 29 compounds with non-significant toxicity or without toxicity as a best ligand. Amongst 29 top ligand the best MM-GBSA dG bind was found to be for FFG7 (docking score = −7.2; MM-GBSA dg bind score = −44.5) followed by FFG2 (docking score = −6.9; MM-GBSA dg bind score = −40.3), which ensures the relative binding energy (Table 1). Furthermore, the output of ROC plot reveals that XP docking performance was better than SP in identifying the actives from decoy.
While the result of docking calculation using Autodock4.2 to find out the binding site of FFG2 and FFG7 along with reference ligand on TR dimer reveals that both these ligands competitively bind at the interface of chain A and chain B of the TR dimer. FFG2 has shown lower binding energy (−12.16 kcal/mol) as compared to FFG7 (−11.42 kcal/mol). While the binding energy for clomipramine (−11.29 kcal/mol) and 4-(4,4,8-Trimethyl-7-oxo-3-oxabicyclo[3.3.1]non-2-yl)-benzoic acid methyl ester (−8.64 kcal/mol) was higher than FFG2 and FFG7. In the binding pocket, both the febrifugine ligand show common hydrogen bonding interactions with Glu-466 and Glu-467. Thus, these interactions with the ligand might be changing the active site structure of TR dimer. The electrostatic potential calculations over TR dimer with febrifugine docked conformation shows that both ligands share the same binding site. We found that the binding energy for FFG2 and FFG7 decreases when the former was performed as compared to Glide ligand docking. The MD simulations of FFG7 and FFG2 showed the overall stability of protein ligand complex. The overall RMSF below 2.5 Å for both FFG7 and FFG2, clearly indicate the stability of complex to be analysed further for in vitro and in vivo parasite inhibition assays.
Therefore, screening and better performance of FFG2 and FFG7 on all above parameter viz. ligand-based molecular docking, ADMET properties, QikProp properties and MM-GBSA. MD reveals that these two compounds may be an effective drug candidate against VL infection.