Molecular Informatics of Trypanothione Reductase of Leishmania major Reveals Novel Chromen-2-One Analogues as Potential Leishmanicides

Trypanothione reductase (TR), a flavoprotein oxidoreductase is an important therapeutic target for leishmaniasis. Ligand-based pharmacophore modelling and molecular docking were used to predict selective inhibitors against TR. Homology modelling was employed to generate a three-dimensional structure of Leishmania major trypanothione reductase ( Lm TR). A pharmacophore model used to screen a natural compound library generated 42 hits, which were docked against the Lm TR protein. Compounds with lower binding energies were evaluated via in silico pharmacological profiling and bioactivity. Four compounds emerged as potential leads comprising Karatavicinol (7-[(2E,6E,10S)-10,11-dihydroxy-3,7,11-trimethyldodeca-2,6-dienoxy]chromen-2-one), Marmin (7-[(E,6R)-6,7-dihydroxy-3,7-dimethyloct-2-enoxy]chromen-2-one), Colladonin (7-[[(4 a S)-6-hydroxy-5,5,8 a -trimethyl-2-methylidene-3,4,4 a ,6,7,8-hexahydro-1 H -naphthalen-1-yl] methoxy]chromen-2-one), and Pectachol (7-[(6-hydroxy-5,5,8 a -trimethyl-2-methylidene-3,4,4 a ,6,7,8-hexahydro-1 H -naphthalen-1-yl)methoxy]-6,8-dimethox-ychromen-2-one) with good binding energies of (cid:1) 9.4, (cid:1) 9.3, 8.8, and (cid:1) 8.5 kcal/mol, respectively. These compounds bound effectively to the FAD domain of the protein with some critical residues including Asp35, Thr51, Lys61, Tyr198, and Asp327. Furthermore, molecular dynamics simulations and molecular mechanics Poisson-Boltzmann surface area (MMPBSA) computations corroborated their strong binding. The compounds were also predicted to possess anti-leishmanial activity. The molecules serves as templates for the design of potential drug candidates and can be evaluated in vitro with optimistic results in producing plausible attenuating infectivity in macrophages.


Introduction
Leishmaniasis is a disease caused by a single-cell eukaryotic parasite of the Leishmania species. This protozoan parasite causes a substantial level of morbidity and mortality. Leishmania has a digenetic life cycle [1]. In mammals, the parasite colonizes macrophages, transforming into intracellular amastigotes. The parasite has an adaptive way to life conditions. The amastigotes tolerate low pH and are hydrolase resistant [2].
Trypanothione is a major product of the trypanothione biosynthesis pathway in trypanosomes which is crucial in maintaining cellular redox potential and is essential for the parasite's survival. This molecule is catalyzed by so many enzymes for which Leishmania major trypanothione reductase (LmTR, E.C. 1.6.4.8) plays a critical role in the biosynthetic pathway. TR reduces trypanothione (T[S] 2 ) to dithiol (T[SH] 2 ). They catalyze the transfer of electrons from NADPH to their specific substrate via an FAD prosthetic group [3]. The reduced form is critical in regulating oxidative stress by reacting with reactive oxygen species (ROS) that are produced by the macrophage. T[SH] 2 is not only needed for detoxification of peroxides but also required for the synthesis of DNA precursors, homeostasis of ascorbate, sequestration and export of thiol conjugate [4].
Trypanothione reductase is a member of the disulphide oxidoreductase family of enzymes. It has an analogue in the human body, glutathione reductase (GR) which also carries out oxidoreductive reactions. But LmTR does not process GSSG and host GR does not reduce T[S] 2 [5,6]. The ascribed reasons for targeting LmTR include the following: (i) trypanothione reductase is a key enzyme in regulating a reducing environment aiding in disease pathogenesis, (ii) this parasite does not depend on the host for reduced trypanothione, (iii) it has a less close known homologous protein in humans; (iv) the availability of template homologs for modeling purpose; and (v) moreover, in vitro trials have proven LmTR to be a good therapeutic target [7].
Computer-aided drug designing is an in-silico approach for drug discovery that combines computational and pharmaceutical research [10]. This application helps in spanning the drug discovery pipeline and helps to speed up and rationalize the drug design process while reducing costs [11]. Ehrlich in 1909 first defined the term pharmacophore as 'a molecular framework that carries (phoros) the essential features responsible for a drug's (pharmacon) biological activity [12]. These features are essential functional groups of atoms in a three-dimensional position that interact with a receptor. Ligand-based drug design can be performed in association with molecular docking. These methods can be combined to identify a number of new hit compounds with potent inhibitory activity and to understand the main interactions at the binding sites. Appropriate use of these methods can improve the ability to identify and optimize hits and confirm their potential to serve as scaffolds for producing new therapeutic agents [13].
Drugs currently used for the treatment of human leishmaniasis are toxic, having severe adverse reactions which limit their use. Aside this includes, increase in resistance by the parasite, high cost of available drugs, lack of efficacy against VL \HIV co-infections with standard chemotherapy, and the development of a single drug or formulation for all forms of leishmaniasis [14][15][16][17]. Therefore, the development of novel, effective drugs with reduced side effects, is still a major priority for health researchers, in spite of many compelling research reports published on antileishmanial agents in the last 10 years [18]. In this study, in silico method of identifying leads was used incorporating the knowledge of pharmacophore and virtual screening to arrive at lead molecules. The overall goal of the study was to predict with high degree potent selective inhibitors of LmTR from the African Natural Product Database (AfroDb) and North African Natural Product Database (NANPDB).

Protein homology modeling
The protein sequence of trypanothione reductase of Leishmania major was retrieved from the NCBI database with the accession number XP_001687512.1, having 491 amino acids [19]. Using the Basic Local Alignment Search Tool (BLAST), the query sequence was compared to known structures which generated structures similar to LmTR. Modeller 9.20 [20] was used for modeling the structure of LmTR.

Active site prediction and quality assessment
To predict the active site, the protein was submitted to CASTp 3.0 [21,22]. The predicted active site was corroborated via a blind docking process using AutoDock Vina within PyRx version 0.9.7 [23,24]. The active region was confirmed by the 'Toggle selection of Spheres' function which highlighted predicted residues from CASTp 3.0. Binding pocket was also viewed with PyMOL v2.0.0 [25]. The quality of the modeled protein was assessed by some quality measure tools. This included PROSA which determines the quality of experimentally solved structures and theoretical models in protein engineering by comparing these to that of experimentally solved protein structures in the PDB database [26]. Verify3D was used to validate the three-dimensional structure of the model [27]. PROCHECK, a quality assessment tool was also used to check the stereochemical properties of the model by generating a Ramachandran plot [28]. ProQ was also used to carry out further validation. ProQ predicted protein quality based on the LGscore and MaxSub scores [29].

Energy minimization of protein target
The modeled LmTR was energy minimized using GROMACS 5.1.1 [30,31]. Simulations were performed with the force field, GROMOS96 43a1. The system was solvated using an equilibrated SPCE216 water model. The charged protein had a net charge of À9 which was neutralized with an equal amount of Na + ions. Energy minimization of the protein was then carried out.

Pharmacophore generation
Pharmacophore model, virtual screening, and molecular docking studies were performed to find novel LmTR inhibitors. The ligand-based structural design incorporates the absence of the macromolecular structure by generating pharmacophore models from a set of ligands. This method takes advantage of the conformational flexibility of the ligands [32]. The active compounds, CNQB and PNTPC were retrieved with IDs CID1323435 and CHEMBL242165 from PubChem (https://pubchem.ncbi.nlm.nih.gov)and CHEMBL (https://www.ebi.ac.uk/chembl) [33,34], respectively and were used to generate the pharmacophore. These ligands served as a training set for pharmacophore generation. All customized settings were kept in default.

Library preparation and pharmacophore validation
LigandScout 4.3 was used for screening a total of 5813 compounds including 885 AfroDb entries found in ZINC database [35] and 4928 NANPDB compounds [36]. The two actives were converted to SMILES format and submitted to the Directory of useful decoys and enhanced (DUD-E) database to generate decoys for the screening [37]. A total of 100 decoys were generated and used as a decoy library. The libraries were then converted into a .ldb file format. The reliability of the pharmacophore model was validated by the area under the receiver operating characteristic (ROC) curve (AUC) [38] using two descriptors, selectivity and sensitivity.

Screening with pharmacophore
LigandScout 4.3 allows in-silico screening of compound libraries using pharmacophore models as filter criteria. Database of active ligands was selected and marked in green with decoys marked in red. The screening process was initiated to generate hits corresponding to the pharmacophore model. These compounds were saved in an "sdf" file format to be used in the docking process.

Docking of LmTR
The generated hit compounds were uploaded into PyRx (Version 0.9.6) [23]. The energy of the ligands was minimized using Universal Force Field (UFF) option in Open Babel incorporated in PyRx prior to docking. This was done to obtain 3D ligand structures which constitute atomic elements that have proper bond lengths between their atoms [39]. Ligands were converted to PDBQT format using AutoDock Vina embedded in the PyRx. Predicted active site residues were selected within a grid box of dimensions X: 42.39 Å, Y: 35.47 Å, and Z: 31.05.14 Å; and centre X: 28.58 Å, Y: 57.09 Å, and Z: À2.24 Å within the AutoDock Vina environment of PyRx for docking process.

Docking validation with AUC
Validation of the algorithm used for the docking process was carried out by generating an AUC plot. Decoys of five known inhibitors in complex with Leishmania trypanothione reductase of several species of Leishmania which included 2,3,4,6-tetra-O-acetyl-1-thio-beta-beta-D-glucopyranose (Auranofin); 4- were retrieved from DUD-E to generate an AUC curve. The result was correlated and plotted using their respective binding energies as the only variable via easyROC (Ver. 1.3) [40].

Docking validation via superimposition and alignment
Superimposition of the crystallographic ligand and re-docking poses was used as a means of validating docking. The five crystallographic ligands in complex with Leishmania trypanothione reductase were removed from the co-crystallized complex and re-docked. The pdb files of the re-docked complex were uploaded into PyMOL together with the solved complex from the Protein Data Bank. LigAlign [41] was then used to calculate the root mean square deviation between the superimposed re-docked and co-crystalized ligands. Superimposition also allowed the identification of critical overlapping residues via LigPlot + . The FAD molecule from the template selected for modeling (PDB ID: 2JK6) was also extracted and docked in LmTR. A comparative study of the original FAD-2JK6 complex and FAD-LmTR complex was done using LigPlot + .

Identification of lead compounds
To identify lead compounds the binding energy, molecular bond interactions, pharmacological, and physiochemical properties were considered. This step helped to filter generalized hit compounds. These compounds in SMILES format were submitted to SwissADME [42], which calculates the corresponding ADME (absorption, distribution, metabolism, and excretion) properties of the compounds. Hydrogen bond interactions of the ligand-protein complexes were studied using LigPlot + and PyMOL.
The hit compounds were physiochemically profiled to identify their druglikeness and solubility in water. Lipinski's rule of 5 was used as a metric to narrow down druggable compounds [43]. Pharmacokinetic properties of predicted compounds were determined in silico. This included cytochrome inhibition, P-glycoprotein (P-gp) substrates, gastrointestinal (GI) absorption, and the bloodbrain barrier (BBB) permeant.

Prediction of activity spectra for substances (PASS) for leads
PASS assesses the probability that a compound has a suspected biological activity [44]. It has been well known that each substance has a wide spectrum of biological activities as evident from some new uses of many old drugs. The SMILES format of the leads were submitted to Way2Drug.com [45] to predict possible biological activity.

Molecular dynamics simulation and MM-PBSA calculation of protein-ligand complexes
By employing GROMACS 2018 [31], the chain A of LmTR and the LmTR-ligand complexes were subjected to molecular dynamics simulations. Ligand topologies were generated via PRODRG which were converted to .gro files. Solvation of each of the protein-ligand complexes in a dodecahedron box was followed by neutralization of the output with sodium and chlorine ions. Each complex was minimized using the steepest descent algorithm coupled with the GROMOS43A force field. Equilibration protocol was carried out to restrain and relax protein-ligand positions. The MD production run was carried out for 100 ns. The output file was used in downstream processes to generate root mean-square deviation (RMSD), root mean-square fluctuation (RMSF) and radius of gyration (Rg) plots with Xmgrace (https://plasma-gate.weizmann.ac.il/Grace/). Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) was employed in calculating the free energies of the complexes. MM-PBSA was carried out using g_mmpbsa, which calculates binding energy components and the per-residue energy decomposition [46]. Graphs were generated using R-programming showing energy interactions.

Homology modeling
L. major trypanothione reductase was modeled with an appreciable degree of accuracy and validated via several bioinformatics tools ( Figure 1A). The protein sequence of L. major with the accession number XP_001687512.1, having 491 amino acids when searched amongst protein homologs in NCBI resulted in a list of similar structures with PDB codes 2JK6, 2X50, 6ER5, 1FEA from Leishmania infantum with their respective percentage identity of 95.72, 95.71, 95.70, and 78.78%, respectively. The similarity of structure and sequence between LmTR protein sequence and the template 2JK6 was 95.72%. This favored its selection for the modeling process. The best model selection was based on the least discrete optimization potential energy (DOPE) score which informs about the energy of the protein.

Active site prediction
In predicting the active site of the protein, results obtained showed 73 binding pockets. Several pockets were identified but the pocket with the largest volume and surface area of 595.278 Å 3 and 924.887 Å 2 , respectively was selected. Larger pockets favor conformational rotation during virtual screening. A total of 80 amino acid active site residues were predicted from CASTp 3.0 ( Table A1). The predicted binding site was visualized using PyMOL ( Figure 1B). The predicted active site was finally confirmed to be FAD binding site. TR has been studied to be a homodimer protein constituting FAD-binding, NADPH-binding, central, and interface domains [5]. The predicted site was concluded to be an FAD site in LmTR by comparative studies between the FAD-2JK6 and FAD-LmTR complexes ( Figure A1-A). The study showed common hydrogen bonding residues including Ser14, Gly15, Arg287, and Thr335 interacting with FAD. These residues correlate with that which was predicted by CASTp 3.0 including other several common residues participating in hydrophobic interactions such as Gly13, Asp35, Ala46, Gly50, Thr51, and Cys52 ( Figure A1-A).

Structural validation and quality prediction
PROSA was used to determine the quality of the structure by comparing it to protein structures that are experimentally solved in the PDB database. It validated the model based on the "quality score or z-score" with a value of À11. 68. The z-score shows whether the predicted model is of X-ray or NMR quality with regards to the amino acid residues length. This z-score value showed that the modeled protein fell in the range of proteins solved experimentally by X-ray crystallography ( Figure A2-A). Verify3D validated the 3D structure of the model. A good 3D structure is expected to have at least 80% of its amino acids to have scored greater than or equal to 0.2 in a 3D/1D profile. This model passed with an appreciable result of 91.65% of residues having a score greater or equal to 0.2 ( Figure A2-B). Further validation with PROCHECK resulted in the generation of a Ramachandran plot ( Figure A3). The plot described the rotations of the polypeptide backbone around the bonds between N-Cα (Phi, φ) and Cα-C (Psi, ψ). The plot allowed the viewing of the distribution of these torsional angles taking into consideration the allowed rotations and rotations that are unfavored which can result in a collision or steric hindrance. A protein with 90% of its residues in the most favorable region is considered a good model. ProQ predicted the quality of protein based on the MaxSub and LGscore scores [29]. The predicted LG score was 6.447 and MaxSub score was 0.520.
LGscore >4 implied that the model was extremely good. Also, a MaxSub score > 0.5 implied that the model was very good [47].

Pharmacophore modeling
The active ligands used as training sets to develop a pharmacophore allowed features similar to the two compounds to be identified and combined into a single geometric function as the basis for the generation of the pharmacophore (Figure 2). Pharmacophore generated utilized features that contributed regions of hydrophobicity and hydrogen bond acceptors incorporated in the model for selective screening. The oxygen from nitrogen dioxide contributed to hydrogen bond acceptors with aromatic and alkene groups contributing to the hydrophobic region ( Figure A4). A number of 10 hypotheses were developed for the model. The best hypothesis with a similarity of 58.12% had a score of 0.8537 and was selected based on the AUC score of 0.99 generated by LigandScout (Figure A5-A). The AUC score was used as a metric to validate how best the pharmacophore model created could distinguish rightly between active compounds and decoys. This intends to reduce false positives and negatives during the screening process.
The screening process was successfully completed with the model which identified 42 compounds that matched the pharmacophore model with a pharmacophore fit score ranging from 55.32 to 57.98 (Table A2).

Docking validation
The docking system was validated using the five inhibitors and 250 decoys generated with DUD-E and further used their binding energies to plot an AUC with a value of 0.702 ( Figure A5-B). AUC value of 1.0 verifies that the prediction of hits obtained from the hypothesis is perfect whereas values of 0.5 and less than 0.70 imply average and moderate random selection respectively [48].
Furthermore, the validation of the molecular docking was also undertaken by aligning the re-docked ligands with their respective co-crystallized complexes taken from PDB. The RMSD values of the alignment between the re-docked ligands and the co-crystalized ligands in complexes of 2YAU, 4APN, 5EBK, 6ER5, 6I7N were 1.483, 3.020, 1.920, 2.712, and 2.465 Å, respectively. Only two of the RMSD values (5EBK and 2YAU) of the alignments were below 2.0 Å, which is considered the threshold for good alignment.
Superimposition also validated the accuracy of docking at the predicted active site. The FAD molecule from the template selected for modeling, 2JK6 when extracted and docked in LmTR, showed similar surrounding residues in the pocket of FAD-2JK6 and FAD-LmTR complex. Ligand alignment of these two complexes gave an RMSD of 3.291 Å which is above the expected threshold. Nonetheless, the FAD-LmTR docking simulated a pose that showed common residues such as Ser14, Gly15, Arg287, and Thr335 taking part in hydrogen bonding in these complexes ( Figure A1). All these verified that the docking system performed very well in docking ligands to the active site.

Virtual screening of pharmacophore hits
When docking validation was verified, molecular docking was carried out. Molecular docking predicted various conformations of each ligand in the binding site of the LmTR. Compounds were selected based on their binding energies. The binding energies gave a theoretical value that relates the affinity of the ligand to the protein model. The results of the respective pharmacophore fit scores and binding energies are well documented ( Table A2). The 42 compounds obtained were Leishmania -One of the Most Diverse Anthroponotic Parasitic Diseases narrowed down to 11 by considering their binding energies and pharmacophore fit scores ( Table 1). With respect to these studies to find inhibitors of the FAD binding site for which FAD is a prosthetic group that already binds tightly to the catalytic site, it was reasonable to select the compounds with binding energies below À9.0 kcal/mol or relatively closer to À9.0 kcal/mol which can provide a plausible binding advantage when acting as inhibitors at the site. The compound ZINC95486081 had the lowest binding energy of À9.8 kcal/mol and the highest was observed from evoxine as À6.2 kcal/mol (Table A2)

Protein-ligand interaction
Molecular interaction studies are important for understanding the mechanism of biological regulation at the molecular level and as such also provides a theoretical basis for drug design and discovery [49,50]. Hydrogen and hydrophobic interactions are key players in stabilizing energetically favored ligands, in an open conformational environment of protein structures [29]. The intermolecular interaction and bond lengths of these 11 compounds were observed. The compound which showed the highest binding affinity, ZINC95486081 formed two hydrogen bonds with residues Lys60 and Ser178 with respective bond lengths of 2.92 Å and 2.87 Å. Five compounds including MTPA, Karatavicinol, Marmin, Betaxanthin and ZINC38658035 had four residues as the highest number of residues partaking in hydrogen bonding. MTPA formed hydrogen bonds with residues Thr51, Thr293, Asp327, and Ser14. Karatavicinol on the other hand formed hydrogen bonds with Thr51, Arg287, Thr293, and Asp327 (Figure A6-A). Marmin formed hydrogen bonds with Val34, Thr51, Thr160, and Thr335 (Figure A6-B) exhibited by Marmin with Thr160. Betaxanthin showed the highest pharmacophore fit score of 57.51 followed by Pectachol (57.18). 13-Hydroxyfeselol showed the lowest fit score of 55.62.

Pharmacological profiling
To identify lead compounds, the binding energy, molecular bond interactions, pharmacological, and physiochemical properties were considered. This step helped to filter generalized hit compounds. The top 11 compounds were profiled in silico to characterize compounds with drug-likeness and good water solubility. Lipinski's rule of 5 was used as a metric to narrow down druggable compounds. The rule factors in the compound's molecular weight which should not exceed 500 g/mol, hydrogen bond donors must not be more than 5, log-p value must be less than or equal to 5 and hydrogen bond acceptors must not be more than 10. All the 11 top hits passed Lipinski's rule of 5 with a good bioavailability score of 0.55 (Table A3).

Pharmacokinetic properties
Further filtering analysis subjected all 11 pharmacophore hits to pharmacokinetics profiling taking into consideration parameters such as gastrointestinal (GI) absorption, blood-brain barrier (BBB) permeation, the permeability of glycoprotein (Pgp), and cytochrome P450 (CYP). Physical parameters such as drug solubility may affect oral bioavailability but in most cases, the major determining factors are likely to be metabolism by CYP and absorption at the intestinal level [51]. CYP3A4 has been known to be responsible for the metabolism of about 50% of all drugs [52] and therefore inhibition of cytochrome can affect oxidation of substrates in cells. Absorption of drugs in the intestine if found high favors the efficacy of the compound as a drug. Multi-drug resistance transporters, such as P-glycoproteins, are essential for many cellular processes that require the transport of substrates across cell membranes [53]. Compounds that are P-gp substrates may face continual efflux which can affect the efficacy of drugs. The blood-brain barrier (BBB) prevents the brain uptake of most pharmaceuticals [54]. This is a disadvantage to neurological diseases but would be of merit since the disease of study is not related to the brain. Compounds that cross the blood-brain barrier may elucidate unwanted biological activities that could be dangerous to health. Therefore, the negative inference would be good for the compound. ZINC95486081 was predicted to show inhibition to three CYP isoenzymes. Karatavicinol, ZINC38658035, and Marmin excelled with an appreciable result (Table A4). For the purpose of narrowing down leads with potential for further computational analysis, compounds with low gastrointestinal absorption were side-lined. This included Taccalin and Betaxanthin.

Prediction activity spectra for substance (PASS)
The biological activity of the selected drug-like candidates was then evaluated using PASS. It is well known that each substance has a wide spectrum of biological activities as evident from some new uses of many old drugs. This allows the tool to utilize this information to predict biological activities based on their probable activity (Pa) and probable inactivity (Pi). When Pa is greater than Pi (Pa > Pi), the compound is likely to possess the predicted biological activity [55,56]. PASS predicted Karatavicinol, Marmin, Colladonin and Pectachol to be potential antileishmanial agents (Table A5) Pectachol with a Pa of 0.694 and Pi of 0.009. Betaxanthin had no prediction as an antileishmanial agent.

Selection of lead compounds for MD and MM-PBSA analysis
The various lead compounds were considered for selection based on the criteria above. ZINC95486081 and MTPA compound although had high binding energy trailed in pharmacokinetic properties and showed Pa less than 0.500. We eliminated Taccalin and Betaxanthin because of their low GI absorption and low Pa values (Tables A3 and A4). Compounds predicted with good probable activity for antileishmanial activities included Karatavicinol, Taccalin, Marmin, 13hydroxyfeselol, Colladonin, Feselol and Pectachol. A literature search revealed Feselol to have antiprotozoal activity against Trypanosoma brucei (IC 50 [57]. That of Feselol and hydroxyfeselol was eliminated because it had been worked on experimentally and also extracts from Ferula genus which is known to exhibit antiviral, antibacterial, and antileishmanial properties [58]. Marmin and Pectachol presented optimistic potential to look into. The high number of hydrogen bonds and binding energy allowed for the inclusion of Karatavicinol in the MM-PBSA to observe their stabilization with this protein target. Also, the high prediction of Colladonin as a compound with probable activity also increased the chance of this compound for MM-PBSA analysis. Therefore, Marmin Pectachol, Karatavicinol, and Colladonin were considered for further analysis in molecular dynamics.

Molecular dynamic simulation of protein-ligand complex
Molecular dynamics simulation allowed the early view of proteins as relatively rigid structures to be replaced by a dynamic model in which the internal motions and resulting conformational changes play an essential role in its function [59]. An RMSD plot generated after molecular dynamics simulation showed a deviation of about 0.25 Å (Figure A7). Further scrutinized with molecular dynamics simulations gave the protein a dynamic dimension to its 3D structural form producing a realistic environment for the ligand interactions that were carried out in the docking process.
Molecular dynamics simulations can also capture a wide variety of important biomolecular processes, including conformational change, ligand binding, and protein folding [60]. The stability of docked protein-ligand complexes was determined by their (RMSD) plots generated from the MD simulation output file. The backbones of the four complexes were observed to be stable over time (Figure 3). The fluctuations of the protein-ligand complexes were analyzed within the system to check for movement and structural stability during the course of the simulation. These movements and stability are significant for the complex functioning inside living systems. The backbone of the LmTR-ZINC8782981 complex showed the greatest stability with an average RMSD of 0.25 nm amongst all the complexes. The LmTR-CHEMBL1277380 complex was fairly stable with RMSD of 0.4 nm. Amongst the four leads, Pectachol was found with the lowest RMSD of 0.3 nm. In terms of stability, the compound Marmin and Colladonin proved to be very much stable around 0.73 nm over the production time of 100 ns. The RMSD of Karatavicinol was observed to show stability from 0 to 60 ns and increased its RMSD to 0.5 nm from 60 to 100 ns.
The flexibility of residues contribution by the LmTR was assessed by the root mean square fluctuation (RMSF). RMSF indicates the flexibility of different regions of a protein, which can be related to crystallographic B factors [61]. The results of the RMSF plots showed consistency for the docked complexes ( Figure A8). The highest fluctuations exhibited was observed around residue numbers 70-90, with Karatavicinol and CHEMBL1277380 showing higher fluctuation levels followed by ZINC8782981 and Marmin complexes. Other regions where good fluctuations were observed include residues between 395 and 410 and 465-480. Overall, Marmin showed more fluctuations around most residues with a distinct difference at residue numbers 260 and 310.
The compactness of the complexes over simulation time is determined by the Rg. If a protein is folded well, it will likely maintain a relatively steady value of Rg, whereas its value will change over time if the protein unfolds [62]. Rg values of all complexes indicated stable complexes over 100 ns (Figure A9). The Rg graph showed most compounds experienced a fairly stable Rg. Marmin experienced the lowest Rg value around 2.33 nm compared to other complexes. This was followed by Colladonin, Pectachol and Karatavicinol with Rg values of 2.37, 2.42, and 2.45 nm, respectively. Between the known inhibitors, CHEMBL1277380 was observed to have an average Rg value of 2.46 nm whilst ZINC8782981 showed the average highest value of around 2.5 nm. Inferring from the Rg graph, the compactness of the LmTR-Marmin, -Colladonin and -Pectachol complexes were maintained after complex formation.

Other energy terms
Van der Waals forces, electrostatic and polar solvation energies, and SASA are relevant energy terms contributing to the overall free binding energy of the complex. The van der Waals energy refers to the weak attraction existing between the intermolecular forces. The van der Waals energy observed in our study showed Karatavicinol and CHEMBL1277380 to have the lowest and highest energy of À228.565 and À 171.823 kJ/mol, respectively. Colladonin, Marmin, and Pectachol also showed relatively low van der Waals energy of À189.289, À189.229, and À 209.538 kJ/mol, respectively as compared with ZINC8782981 with À222.123 kJ/mol. Electrostatic energy refers to the potential energy of a system consisting of different electric charges [65][66][67]. The lowest electrostatic energy was exhibited by Marmin (À386.401 kJ/mol) followed by Pectachol (À286.260 kJ/mol), and Colladonin (À249.067 kJ/mol). Karatavicinol and the other two inhibitors were observed with high electrostatic energy ( Table 2). Some studies have observed that van der Waals and electrostatic forces contribute favorably to the energetics of binding along with simulations that favor the binding of complexes [66,68].
Polar solvation energy on the other hand refers to the electrostatic interaction that exists between the solute and the continuum solvent [69]. The highest polar solvation energy amongst the leads was exhibited by Marmin (484.074 kJ/mol) and the lowest by Karatavicinol (227.483 kJ/mol). Solvent accessible surface area (SASA) energy was calculated after MD. This represents the non-polar solvation energy [69]. This energy measures the interactions that exist between the complex and the solvents. Amongst the leads, Karatavicinol obtained the lowest SASA energy followed by Pectachol, Colladonin, and Marmin ( Table 2). Relative to these were the low SASA energies of the inhibitors ZINC8782981 and CHEMBL1277380.
The total contribution of these energies enabled the final estimation of the free binding energies in the complexes ( Table 2). The lowest free binding energy contributing to more stability of the protein-ligand complex was observed by Marmin (À109.114 kJ/mol). Next amongst the four complexes was Pectachol (À63.487), followed by Karatavicinol (À57.644 kJ/mol), and Colladonin (À48.936 kJ/mol). The low binding energy of Marmin was much closer to that of CHEMBL1277380 (À111.732 kJ/mol) with that of Pectachol higher than ZINC8782981 (À54.399 kJ/ mol). These energies address the potential of Marmin and Pectachol to bind most effectively at the active site of LmTR. LmTR-Marmin's free binding energy correlated with the low binding energy (À9.3 kcal/mol) from docking. That of Pectachol showed a good free binding energy than that obtained from docking. This was better than that of Karatavicinol and Colladonin ( Table 1).

Exploring possible implications and structure similarities of predicted leads
Karatavicinol and Marmin had lower binding energies of À9.4 and À 9.3 kcal/ mol, respectively, as compared to Colladonin (À8.5 kcal/mol) and Pectachol (À8.5 kcal/mol). These binding energies are closer to that of FAD (À9.0 kcal/mol) for which can possibly compete in binding at the FAD domain. These compounds were concluded to have drug-likeness by satisfying Lipinski's rule of 5. They also do not pass the blood-brain barrier which is good. Also, Marmin and Karatavicinol checked false for p-glycoprotein substrate. This gives the compounds an advantage to maintain their concentrations in cellular level to maximize efficacy. Pectachol and Colladonin however were implicated as P-gp substrates. These predicted preferable properties can favor their lead likeness and chances of going a long way in experimental studies. The four lead compounds were predicted as antileishmanial compounds. The four leads are confirmed not to be already existing antileishmanial drugs by structural similarity searches in www.DrugBank.ca but rather observed to be analogues of chrome 2-one. In regard to this, studies over the years have however shown some novel compounds such as 7-{[(2R*)-3,3-dimethyloxiran-2-yl] methoxy}-8-[(2R*,3R*)-3-isopropenyloxiran-2-yl]-2H-chromen-2-one and 7methoxy-8-(4-methyl-3-furyl)-2H-chromen-2-one against Leishmania donovani with EC 50 of 9.9 and 10.5 μg/mL, respectively [70]. These tested compounds with antileishmanial effect tend to be analogues of chromen-2-one. We emphasize that Karatavicinol is not a unique lead compound since it has already been experimented on other Leishmania species excluding L. major [57]. But the study identified it via these computational processes and therefore would report it as a potential compound against L. major. This augments the fact that the computational drug discovery pipeline has an optimistic potential of yielding good candidates for experimental work. Colladonin on the other hand is an enantiomer of Feselol for which Feselol is experimented as an antileishmanial agent [57]. Marmin also holds a very good potential of being an anti-ulcerative agent [71]. This favors it being a good

24.773
The energy values are presented as mean AE standard deviation kJ/mol. compound for the treatment of cutaneous leishmaniasis. These compounds classified are coumarins and more other studies have reported good antileishmanial activities from this class of compounds [72,73]. This work supports the fact that Karatavicinol, Marmin, Pectachol and Colladonin may possibly exhibit good antileishmanial activity if tested in vitro (Table A6). Further in this study, the interaction of the active site residues with all four lead compounds showed hydrogen bonding with Val34, Thr51, Lys60, Thr160, Ala159, Arg287, Thr293, Asp327, Asn330, Thr335, and Gly376 ( Table 1). Superimposition of the docked 2JK6 and co-crystallized revealed common residues such as Ser14, Gly15, Arg287, and Thr335 ( Figure A1). These residues can be observed to be unique to the FAD domain of LmTR in anchoring the FAD molecule. Comparing these residues to the hydrogen bonding residues from the four leads shows that possible interruption of any of these residues can cause conformational changes which might not favor the selective binding of FAD at its domain. Baiocco et al. in 2009 identified Thr335 of trypanothione reductase at the FAD catalytic site of L. infantum [74]. They proposed that the FAD molecule binds tightly to the protein and orients itself towards the hydride transfer region of the active site by hydrogen bonding with specific residues Lys60, Thr335, and His461. Having observed this, interrupting one of these residues can potentially inhibit the reduction of T[S] 2 by interfering with the hydride transfer. These compounds can potentially covey a competitive mode for binding to Thr355 which can affect the hydride transfer reaction in the active site preventing direct inactivation of trypanothione reductase. Other studies with quinone derivatives also have identified Thr355 and Ser14 as unique to the FAD domain of TR [75,76].

Conclusion
Trypanothione reductase has been a well-investigated target essential for trypanosomatids. Its function in controlling oxidative stress in the parasite provided an opportunity to target the trypanothione biosynthesis pathway. A total of 11 hit compounds identified by pharmacophore modeling and virtual screening were filtered to four potential leads by considering their ADME with their molecular interactions in LmTR. MM-PBSA enabled the individual computation of active site residues that contributed significantly to binding. Efficient selective blockade of LmTR with these four coumarin compounds: Karatavicinol (7-[(2E,6E,10S)-10,11dihydroxy-3,7,11-trimethyldodeca-2,6-dienoxy]chromen-2-one), Marmin   LmTR. This shows 91.65% of its amino acid residues with an average 3D-1D score greater than or equal to 0.2, which is a positive inference to the expected 80% of amino acids with 3D-1D score above or equal to 0.2.

18
Leishmania -One of the Most Diverse Anthroponotic Parasitic Diseases Figure A2. Ramachandran plot of the modeled LmTR protein structure. The percentage of residues in the most favored region (red) was 93.6% which is favorable for the protein's stereochemistry. The percentage of residues in the allowed region (yellow) was 6%. Only 0.2% of protein residues (Phe45) showed probable stereochemical hindrance or collision.

19
Molecular Informatics of Trypanothione Reductase of Leishmania major… DOI: http://dx.doi.org/10.5772/intechopen.100594 Figure A3. A 3D geometry of the generated pharmacophore. The nitrogen on the bicyclic ring of CNQB with the oxygen from the nitro group on its purine ring derivative contributed a hydrogen bond acceptor HBA (red sphere). The oxygen from the nitrogen dioxide group on the conjugated benzene in addition to the nitrogen on the fivemember ring of PNTPC also contributed to HBA. Both had an aromatic ring (blue ring) as a common feature (blue ring in yellow 3D sphere) which contributed to hydrophobic interactions and the alkene feature shared amongst them generated the same hydrophobicity.