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

Abstract

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 (LmTR). A pharmacophore model used to screen a natural compound library generated 42 hits, which were docked against the LmTR 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-[[(4aS)-6-hydroxy-5,5,8a-trimethyl-2-methylidene-3,4,4a,6,7,8-hexahydro-1H-naphthalen-1-yl]methoxy]chromen-2-one), and Pectachol (7-[(6-hydroxy-5,5,8a-trimethyl-2-methylidene-3,4,4a,6,7,8-hexahydro-1H-naphthalen-1-yl)methoxy]-6,8-dimethoxychromen-2-one) with good binding energies of −9.4, −9.3, 8.8, and −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.

1. 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]₂) to dithiol (T[SH]₂). 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]₂ 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]₂ [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].

Several inhibitors have been screened against this enzyme causing a reduction of infectivity and decreased capacity of the parasite to survive within intracellular macrophages. Potent compounds, such as 7-chloro-4-nitro-5-quinazolin-4-ylsulfanyl-2,1,3-benzothiadia-zole (CNQB) and 4-phenyl-5-(4-nitro-cinnamoyl)-1,3,4-thiadiazolium-2-phenylamine-chloride (PNTPC) with IC₅₀ values 0.58 and 1.63 μM, respectively have already been tested in an in vitro assay against trypanothione reductase of trypanosomatids [8, 9].

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).

2. Materials and methods

2.1 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.

2.2 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].

2.3 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.

2.4 Pharmacophore modeling and screening

2.4.1 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.

2.4.2 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.

2.4.3 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.

2.5 Molecular docking

2.5.1 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.

2.5.2 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-[[1-(4-ethylphenyl)-2-methyl-5-(4-methylsulfonylphenyl)pyrrol-3-yl]methyl]thiomorpholine (CHEMBL1277380); 6-sec-butoxy-2-[(3-chlorophenyl)sulfanyl]-4-pyramidinamine (RDS); 2-(diethylamino)ethyl4-((3-(4-nitrophenyl)-3-oxopropyl)amino)benzoate (ZINC8782981); {N}-(4-azanylbutyl)- ∼ {N}-(2-azanyl-2-oxidanylidene-ethyl)-7-(3-azanyl-3-oxidanylidene-propyl)-4-(dimethylamino)-2-(2-naphthal en-2-ylethylamino)pyrrolo[2,3-d]pyrimidine-6-carboxamide (H6H) 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].

2.5.3 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⁺.

2.6 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 drug-likeness 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 blood–brain barrier (BBB) permeant.

2.7 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.

2.8 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.

3. Results and discussion

3.1 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.

Figure 1.

(a) 3D structure of modeled protein in cartoon representation. A monomer of LmTR with alpha helices represented in red, the beta sheets in yellow and the loops in green. (b) Surface representation of the LmTR. Active site (in red) is the FAD binding region.

3.2 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 Å³ and 924.887 Å², 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).

3.3 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].

3.4 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.

Figure 2.

2D structure of PNTPC (left) and CNQB (right) showing shared features. PNTPC and CNQB showed two hydrogen bond acceptors. Both compounds had an aromatic ring as a common feature. These served as active ligand for the training set.

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).

3.5 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.

3.6 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 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).

Predicted ligands	Pharmacophore fit score	Binding energy/(kcal/mol)	Hydrogen bond Residues and length (Å)	Hydrophobic bond interacting residues
ZINC95486081	55.95	−9.8	Lys60 (2.92); Ser178 (2.87)	Gly13, Gly15, Asp32, Ala46, Thr51, Cys52, Val55, Gly56, Ala159, Thr160, Tyr198, Arg287, Asp327, Met333, Leu334, Thr335
MTPA	56.37	−9.4	Thr51 (3.12); Thr293 (2.89); Asp327 (3.13); Ser14 (3.17)	Gly11, Gly13, Gly15, Asp32, Val36, Ala46, Gly50, Cys52, Val55, Gly56, Ala159, Thr160, Tyr198, Arg287, Met333, Leu334, Thr335
Karatavicinol	56.5	−9.4	Thr51 (3.28); Arg287 (3.00); Thr293 (3.01, 2.97); Asp327 (2.82)	Ser14, Gly11, Gly13, Val34, Asp35, Val36, Gly50, Cys52, Ala46, Phe126, Gly127, Ala159, Thr160, Gly161, Tyr198, Arg290, Met333, Ala338
Taccalin	56.42	−9.4	Ser14 (3.21); Ala365 (3.25)	Gly13, Thr51, Cys52 Gly56, Cys57, Lys60, Gly161, Ile199, Thr198, Gly326, Met333, Leu334, Thr335, Ala338
Marmin	56.18	−9.3	Val34 (3.17); Thr51 (2.99, 3.04); Thr160 (2.83); Thr335 (2.87)	Leu10, Gly11, Gly13, Ser14, Asp35, Ala46 Gly50, Cys52, Gly127, Ala159, Gly161, Arg290, Leu294, Ala327, Leu334, Ala338
13-Hydroxyfeselol	55.62	−9.1	Val362 (2.79); Thr374 (2.85); Gly376 (3.20, 3.26)	Lys60, Thr198, Gly229, Phe230, Gly326, Leu334, Cys364, Ala365
Betaxanthin	57.51	−8.9	Ser14 (2.85); Cys52 (3.03); Gly127 (3.03, 3.22, 3.30); Thr335	Gly11, Gly13, Val34, Asp35, Val36, Gly50, Cys52, Ala46, Phe126, Ala159, Thr160, Gly161, Tyr198, Arg290, Asp327 Met333, Ala338
Colladonin	55.90	−8.8	Asn330 (2.85)	Lys60, Gly197, Tyr198, Tyr221, Arg287, Phe230, Leu334, Ala365
Feselol	55.63	−8.8	Asn330 (2.87)	Lys60, Gly197, Tyr198, Tyr221, Arg287, Phe230, Leu334, Ala365
ZINC38658035	55.95	−8.7	Tyr198 (3.31); Val362 (2.92); Tyr374 (3.30); Gly376 (2.86, 3.05)	Ile199, Phe230, Gly286 Arg287, Met333, Leu334, Cys364, Cys375
Pectachol	57.18	−8.5	Lys60 (2.93); Gly376 (3.11)	Tyr198, Gly229, Phe230, Val332, Met333, Leu334, Ala365, Val362, Cys364, Val366, Phe367
FAD molecule and inhibitors from 6ER5 and 4APN
ZINC8782981	_	−7.2	Lys60 (3.12); Arg287 (3.08)	Cys52, Cys57, Tyr198, Gly229, Phe230, Val332, Met333, Leu334, Ala365, Val362, Cys364
CHEMBL1277380	_	−8.2	Lys60 (2.89)	Tyr198, Gly229, Phe230, Val332, Met333, Leu334, Ala365, Val362, Gly376
FAD	_	−9.0	Ser14 (3.11); Gly15 (2.96); Ala159 (2.97); Tyr198 (2.87); Arg287 (2.75); Met333 (2.78); Thr335 (2.88, 3.31)	Gly13, Gly50, Thr51, Cys52, Ser162, Gly197, Gly229, Phe230, Asp327, Leu334, Ala338

Table 1.

The table shows parameters involved in the selection of lead compounds. This included pharmacophore fit score, the binding energy and the number of hydrogen and hydrophobic bond interacting residues.

3.7 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). Betaxanthin also bonded with Ser14, Cys52, Gly127, and Thr335. Finally, ZINC38658035 formed hydrogen bonds with Tyr198, Val362, Tyr374, and Gly376. The compound 13-hydroxyfeselol was the only hit that formed three hydrogen bonds with Val362, Thr374, and Gly376. Taccalin and Pectachol formed only two hydrogen bonds with Ser14 and Ala365. On the other hand, Colladonin and Feselol formed the least hydrogen bond residues with Asn330. The shortest bond length of 2.83 Å was 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.

3.8 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).

3.9 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.

3.10 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). Colladonin showed the highest Pa of 0.768 and Pi of 0.006 followed by Taccalin (Pa of 0.711 and Pi of 0.009) and Pectachol with a Pa of 0.694 and Pi of 0.009. Betaxanthin had no prediction as an antileishmanial agent.

3.11 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, 13-hydroxyfeselol, Colladonin, Feselol and Pectachol. A literature search revealed Feselol to have antiprotozoal activity against Trypanosoma brucei (IC₅₀ 8.1 μM), Trypanosoma cruzi (IC₅₀ 8.6 μM), and Leishmania infantum (IC₅₀ 6.8 μM). Similarly, 10′R-karatavacinol has presented activity against T. brucei (IC₅₀ 32.4 μM), T. cruzi (IC₅₀ 9.4 μM), and L. infantum (IC₅₀ 32.4 μM) [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.

3.12 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.

Figure 3.

RMSD values of the LmTR-ligand complexes of the four leads (Karatavicinol, Marmin, Pectachol, and Colladonin) and the two known inhibitors after 100 ns. The complexes in the graph are color coded.

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.

3.13 Evaluation of leads using MM-PBSA

MM-PBSA was employed to calculate free binding energies by per-residue decomposition of the protein complexes. At a quantitative level, simulation-based methods provide substantially more accurate estimates of ligand binding affinities (free energies) than other computational approaches such as docking [63]. Residues contributing binding free energy greater than 5 kJ/mol or less than −5 kJ/mol are considered critical for binding of a ligand to a protein [64]. MM-PBSA results showed only Asp327 amongst the hydrogen bonding residues of Karatavicinol to contribute a per residue decomposition energy of 13.65 kJ/mol. Amino acid residue Asp35 (21.89 kJ/mol) was observed with such greater contribution (Figure A10). The complex of LmTR–Marmin also showed surrounding hydrophobic residues Asp35 (−8.62 kJ/mol), Ala46 (−7.65 kJ/mol), Arg290 (8.83 kJ/mol), and Glu141 (−5.12 kJ/mol) with their energy decomposition to be greater than 5 kJ/mol and less than −5 kJ/mol. The only hydrogen bonding residue that showed a relevant contribution of energy decomposition was Thr51 (−16.36 kJ/mol) (Figure 4). Hydrophobic amino acid residues Lys61 (−5.16 kJ/mol), Tyr198 (−11.29 kJ/mol), Asp327 (5.56 kJ/mol), and Arg331 (6.38 kJ/mol) showed relevant contribution to the total binding energy of the LmTR–Colladonin complex (Figure A11). Moreover, only hydrogen bonding residue Lys60 (13.32 kJ/mol) in LmTR–Pectachol complex showed to be a critical residue in binding. Other surrounding residues contributed substantially to the per residue energy decomposition in the LmTR–Pectachol complex. This included Lys61 (−10.97 kJ/mol), Arg287 (−6.91 kJ/mol), Asp327 (9.30 kJ/mol), Met333 (−6.72 kJ/mol), Leu334 (−8.91 kJ/mol), Lys361 (6.46 kJ/mol), and Cys364 (−6.07 kJ/mol) (Figure A12). Deducing from the substantial contribution of energy per decomposition of residues, we propose Asp35, Thr51, Lys61, Tyr198, and Asp327 to be critical in intermolecular bonding and stabilization of ligands at the FAD active site.

Figure 4.

MM-PBSA plot of the binding free energy decomposition contribution per residue of LmTR–Marmin complex. Coded red lines represent surrounding active site amino acid residues.

3.14 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].

Compound	Van der Waals energy (kJ/mol)	Electrostatic energy (kJ/mol)	Polar solvation energy (kJ/mol)	SASA energy (kJ/mol)	Binding energy (kJ/mol)
ZINC8782981	−222.123 ± 14.568	−52.495 ± 17.662	245.049 ± 33.654	−24.829 ± 0.962	−54.399 ± 20.084
CHEMBL1277380	−171.823 ± 14.173	−2.926 ± 5.485	82.884 ± 16.706	−19.869 ± 1.202	−111.732 ± 16.514
Karatavicinol	−228.565 ± 12.673	−32.345 ± 21.415	227.483 ± 27.305	−24.217 ± 1.100	−57.644 ± 24.019
Marmin	−189.289 ± 16.726	−386.401 ± 30.540	484.074 ± 28.991	−17.498 ± 1.050	−109.114 ± 23.461
Pectachol	−189.229 ± 18.203	−286.260 ± 49.152	430.604 ± 71.136	−18.602 ± 1.308	−63.487 ± 33.289`
Colladonin	−209.538 ± 18.908	−249.067 ± 40.851	427.216 ± 49.348	−17.548 ± 1.122	−48.936 ± 24.773

Table 2.

The energy terms obtained after MM-PBSA analysis of the protein–ligand complexes.

The energy values are presented as mean ± standard deviation kJ/mol.

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).

3.15 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 7-methoxy-8-(4-methyl-3-furyl)-2H-chromen-2-one against Leishmania donovani with EC₅₀ 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 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]₂ 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].

4. 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,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), Pectachol (7-[(6-hydroxy-5,5,8a-trimethyl-2-methylidene-3,4,4a,6,7,8-hexahydro-1H-naphthalen-1-yl)methoxy]-6,8-dimethoxychromen-2-one), and Colladonin (7-[[(4aS)-6-hydroxy-5,5,8a-trimethyl-2-methylidene-3,4,4a,6,7,8-hexahydro-1H-naphthalen-1-yl]methoxy]chromen-2-one) hold the potential to compromise the redox defenses of the parasites by inhibiting the FAD binding region and correspondingly increasing their sensitivity to redox-damage when carried out in in vitro and in vivo studies. Residues such as Asp35, Thr51, Lys61, Tyr198, and Asp327 are suspected to have critical role in the anchoring of FAD which contributes to the formation of reduced T[SH]₂ in the reducing environment of amastigotes.

Acknowledgments

The authors are grateful to the West African Centre for Cell Biology of Infectious Pathogens (WACCBIP) at the University of Ghana for making Zuputo, a Dell EMC high-performance computing cluster, available for this study.

Conflict of interest

The authors declare no conflict of interest.

Leu10, Gly11, Gly13, Ser14, Gly15, Gly16, Val34, Asp35, Val36, Phe44, Ala46, Ala47, Gly50, Thr51, Cys52, Val55, Gly56, Cys57, Lys60, Lys61, Gly125, Phe126, Gly127, Ala128, Arg138, Ser140, Glu141, Pro143, Ala159, Thr160, Gly161, Ser162, Trp163, Pro164, Thr165, Thr177, Ser178, Asn179, Phe182, Tyr198, Ile199, Glu202, Phe203, Met282, Leu283, Ala284, Ile285, Gly286, Arg287, Arg290, Thr293, Leu294, Gln295, Ile325, Gly326, Asp327, Val332, Met333, Leu334, Thr335, Pro336, Val337, Ala338, Ile339, Asn340, Arg355, Thr357, Asp358, His359, Thr360, Lys361,Val362, Ala363, Cys364, Ala365, Phe367, Pro435, Glu436, Ile438, Gln439, Gly442, Ile443, Lys446

Table A1.

Predicted amino acid residues.

Figure A1.

(a) The Z-score of the modeled LmTR (represented in dot) estimated to be −11.68 which is within the range of experimentally determined proteins by X ray method. (b) Verify 3D plot of the modeled protein structure, 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.

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.

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 five-member 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.

Figure A4.

(A) AUC score of 0.99 for the pharmacophore model. Determined at 1, 5, 10, and 100% of the selected database were the AUC and EF values as shown. The median is shown by dotted lines. If the curve is closer to the median it would suggest poor model. (B) AUC score of 0.702 generated for validating the docking system used. It verified the correlation between virtual screening performance and binding site descriptors of protein targets model (LmTR).

Figure A5.

(a) 2D schematic diagram of co-crystallized FAD (PDB ID: 2JK6) and FAD docked in LmTR superimposed together. Similar hydrogen bonding residues include Ser14, Gly15, Arg287, and Thr335. Similar hydrophobic residues in addition confirm the predicted active site. (b) Ligand alignment of co-crystalized FAD and FAD docked in LmTR.

Name	P-fit score	Binding energy (kcal/mol)	Active/decoy	Source database
ZINC95486081	55.95	−9.8	Active	AfroDB
(S)-alpha-methoxy-alpha-trifluoromethyl-alpha-phenylacetate (MTPA)	56.37	−9.4	Active	NANPDB
Karatavicinol	56.5	−9.4	Active	NANPDB
Taccalin	56.42	−9.4	Active	NANPDB
Marmin	56.18	−9.3	Active	NANPDB
3-Hydroxyfeselol	55.62	−9.1	Active	NANPDB
ZINC95486257	55.9	−9.0	Active	AfroDB
Betaxanthin	56.97	−8.9	Active	NANPDB
Coladonin	56.58	−8.8	Active	NANPDB
Feselol	56.41	−8.8	Active	NANPDB
ZINC38658035	55.9	−8.7	Active	AfroDB
Pectachol	57.18	−8.5	Active	NANPDB
ZINC85967928	55.85	−8.4	Active	AfroDB
Polyanthin	56.39	−8.4	Active	NANPDB
ZINC95486047	57.98	−8.3	Active	AfroDB
4′-Methyl gossypetin	56.17	−8.2	Active	NANPDB
2-(nonan-8-one)-4-methoxy-quinoline	56.41	−8.2	Active	NANPDB
Orientin	55.53	−8.1	Active	NANPDB
Kaempferol-3,6-dimethylether-7-glucoside	57.15	−7.8	Active	NANPDB
ZINC95486129	56.43	−7.8	Active	AfroDB
Ethuliaconyzophenone	56.9	−7.7	Active	NANPDB
ZINC95486209	56.55	−7.5	Active	AfroDB
(+)-1,2-bis-(4-hydroxy-3-methoxyphenyl)-propane-1,3-diol [erythro form]	55.9	−7.4	Active	NANPDB
4-Hydroxy-2′,4′-dimethoxy-dihydrochalcone	55.58	−7.4	Active	NANPDB
Drimartol A	56.31	−7.4	Active	NANPDB
Isoarnottinin-4′-O-beta-d-glucoside	55.71	−7.4	Active	NANPDB
4-Beta-hydroxy-6alpha-(4-hydroxy-3-methoxybenzoyl)-7-daucen-9-one	55.93	−7.4	Active	NANPDB
ZINC14686464	56.55	−7.4	Active	AfroDB
6-(3′,4′-dimethoxybenzoyl)-jaeschkeanadiol	57.17	−7.3	Active	NANPDB
ZINC14887523	56.88	−7.3	Active	AfroDB
Orientin-7-methoxide	56.26	−7.2	Active	NANPDB
ZINC14444870	56.35	−7.2	Active	AfroDB
ZINC14689062	56.5	−7.2	Active	AfroDB
1-Dehydrogingerdione	56.05	−7.1	Active	NANPDB
Onopordin	56.27	−7.1	Active	NANPDB
ZINC95486194	56.79	−7.1	Active	AfroDB
Methyl5-(3-4-dihydroxyphenyl)-3-hydroxypenta-2,4-dienoate	55.32	−7	Active	NANPDB
Corniculatusin	56.23	−7	Active	NANPDB
3-(10-acetoxygeranyl)-4-acetoxy-p-coumaric acid	56.14	−7	Active	NANPDB
ZINC00035526	56.66	−7	Active	AfroDB
ZINC00608186	57.08	−6.8	Active	AfroDB
Evoxine	57.32	−6.2	Active	NANPDB

Table A2.

The 42 hits obtained from pharmacophore screening with their respective pharmacophore fit score, binding energies, and data sources.

Figure A6.

2D schematic diagram showing protein-ligands interaction of some leads at the active site of LmTR. (A) LmTR–Karatavicinol interaction profile, (B) LmTR–Marmin interaction profile, (C) LmTR–Pectachol interaction profile, and (D) LmTR–Colladonin interaction profile.

Compound ZINC ID/name	Number of Lipinski’s rules violated	MW (g/mol)	No. HA	No. HD	xLogP	Water solubility (mg/mL)	Log S	Bio. Sc
ZINC95486081	0	382.45	5	2	4.52	Moderately soluble	−5.84	0.55
MTPA	0	470.52	8	0	6.35	Moderately soluble	−6.00	0.55
Karatavicinol	0	400.51	5	2	4.66	Moderately soluble	−4.85	0.55
Taccalin	0	418.48	9	6	−1.45	Moderately soluble	1.66	0.55
Marmin	0	332.39	5	2	2.81	Soluble	−3.52	0.55
13-Hydroxyfeselol	0	400.51	5	2	1.45	Moderately soluble	−5.93	0.55
Betaxanthin	0	370.44	8	7	−1.17	Moderately soluble	−2.11	0.55
Colladonin	0	384.51	4	1	5.76	Poorly soluble	−6.50	0.55
Feselol	0	384.51	4	1	5.76	Poorly soluble	−6.5	0.55
ZINC38658035	0	464.63	6	3	−4.47	Soluble	−3.28	0.55
Pectachol	0	444.56	6	1	5.70	Poorly soluble	−6.70	0.55

Table A3.

Physicochemical profiling of the 11 hit compounds.

All the hits showed good druglikeness. MW, molecular weight; No. HD, number of H-bond donors; Bio Sc, bioavailability score, No. HA, number of H-bond acceptors.

Compound ZINC ID	GI absorption	BBB permeant	P-gp substrate	CYP1A2 inhibitor	CYP2C19 inhibitor	CYP2C9 inhibitor	CYP2D6 inhibitor	CYP3A4 inhibitor
ZINC95486081	High	Yes	Yes	No	No	Yes	Yes	Yes
MTPA	High	No	No	No	No	No	Yes	Yes
Karatavicinol	High	No	No	No	No	No	No	Yes
Taccalin	Low	No	Yes	No	No	No	No	No
Marmin	High	No	No	No	No	No	No	No
13-Hydroxyfeselol	High	No	Yes	No	No	No	Yes	Yes
Betaxanthin	Low	No	No	No	No	No	No	No
Colladonin	High	Yes	No	No	No	No	Yes	No
Feselol	High	Yes	No	No	No	No	Yes	No
ZINC38658035	High	No	Yes	No	No	No	No	No
Pectachol	High	No	Yes	No	No	No	Yes	No

Table A4.

Pharmacological profiling of top 11 compounds characterized by gastrointestinal (GI) absorption, blood brain barrier (BBB) permeant, p-gp substrates, and cytochrome P450 inhibitors.

Four compounds out of the 11 showed an appreciable pharmacological property. This included Karatavicinol, Marmin, Colladonin, and Pectachol.

Lead compounds	Antileishmanial predicted activity
	Pa	Pi
ZINC95486081	0.224	0.168
MTPA	0.263	0.130
Karatavicinol	0.513	0.021
Taccalin	0.711	0.009
Marmin	0.557	0.024
13-Hydroxyfeselol	0.658	0.030
Betaxanthin	—	—
Colladonin	0.768	0.006
Feselol	0.768	0.006
ZINC38658035	0.345	0.074
Pectachol	0.694	0.009

Table A5.

This table shows the 10 top hit compounds and their predicted antileishmanial activity.

Karatavicinol, Taccalin, Marmin, 13-Hydroxyfeselol, Colladonin, Feselol, and Pectachol had a greater positive prediction above 0.5. If 0.5 < Pa < 0.7, the substance is likely to exhibit activity in experiment, but the probability of being a known pharmaceutical agent is less.

Figure A7.

Root mean square fluctuations of six complexes. The complexes are color coded in the graph. Karatavicinol and CHEMBL1277380 experienced the highest fluctuation at around residue number 80. Remaining complexes had similar patterns of fluctuations.

Figure A8.

The radius of gyration (Rg) plots of seven complexes within 100 ns simulation time. The complexes are represented in color code in the graph. Marmin showed the most preferentially well folded protein complex with Rg value of 2.33 nm.

Figure A9.

MM-PBSA plot of the binding free energy decomposition contribution per residue of LmTR–Karatavicinol complex. Coded red lines represent surrounding active site amino acid residues.

Figure A10.

MM-PBSA plot of the binding free energy decomposition contribution per residue of LmTR–Pectachol complex. Coded red lines represent surrounding active site amino acid residues.

Figure A11.

MM-PBSA plot of the binding free energy decomposition contribution per residue of LmTR–Colladonin complex. Coded red lines represent surrounding active site amino acid residues.

Figure A12.

Shows a graph of RMSD of the backbone of atoms in nm versus time in nanoseconds (ns). This graph is a representation of the average distance of the atoms of the residues at the backbone of the target protein. RMSD of 0.25 Å showed deviation from protein backbone.

Table A6.

Structure and IUPAC names of the three novel lead compounds.

A.1 List of abbreviations

ADMET

absorption, distribution, metabolism, excretion and toxicity

AUC

area under curve

CYP

cytochromes P₄₅₀

DUD-E

directory of useful (docking) decoys-enhanced

GROMACS

GROningen MAchine for Chemical Simulations

HPC

high performance computing

ID

identification

Log P

logarithm of the octan-1-ol/water partition coefficient

MD

molecular dynamics

MM-PBSA

molecular mechanics Poisson Boltzmann surface area

Mw

molecular weight

P-gp

permeability glycoprotein

PASS

prediction of activity spectra for substance

PDB

Protein Data Bank

Rg

radius of gyration

RMSD

root mean square deviation

RMSF

root mean square fluctuation

ROC

receiver operating characteristic

SDF

structure data file

SMILES

simplified molecular input line entry system

UFF

universal force field