1. INTRODUCTION
Visceral leishmaniasis, also known as kala-azar or black fever, is a neglected tropical disease caused by the intramacrophage parasite, Leishmania donovani (L. donovani) [1]. This parasite has a digenetic life cycle; the promastigotes form occurs in female sandflies (Phlebotomus sp.) and the amastigotes form grows and multiplies in the macrophages of humans and mammals hosts [2,3]. It primarily affects the internal organs, most notably the bone marrow, liver, and spleen [4]. In 2020, ten countries (India, Kenya, Eritrea, Yemen, China, Brazil, Somalia, Ethiopia, Sudan, and South Sudan) reported more than 90% of global visceral leishmaniasis cases [5]. Given the absence of a viable vaccine, current treatment is limited to a few expensive drugs such as pentavalent antimonials, miltefosine, pentamidine, and amphotericin B [6]. The use of these drugs is also constrained by severe adverse effects, lengthy treatment duration, and parasitic drug resistance [4,6,7]. The current scenario requires the development of new and secure medications to supplement the currently available therapies, thereby compelling the need for this study. Leishmania donovani adenosylmethionine decarboxylase (LdAdoMetDC) of the polyamine pathway has been reported as a potential target for antileishmanial therapy [8-10]. It is an obligatory enzyme present in higher eukaryotes as well as in eukaryotic protozoa trypanosomatids such as Trypanosoma brucei, L. donovani, and other trypanosomatids [11]. Adenosylmethionine decarboxylase (AdoMetDC) is responsible for the irreversible decarboxylation of S-adenosylmethionine, leading to the production of S-adenosyl-5´-(3-methylthio propylamine). This compound, in conjunction with putrescine, acts as a substrate for the enzyme spermidine synthase [12,13]. Spermidine is required for the parasite’s viability, growth, and infectious mammalian stage [8,14,15]. In 2002, the gene encoding AdoMetDC had been cloned and characterized from L. donovani and Leishmania infantum [8,9,16]. Gene deletion studies in mice at very early embryonic stages of L. donovani established that AdoMetDC is an essential enzyme of the polyamine pathway [9,17]. AdoMetDC exhibited high homology with the various trypanosomatid species (62–85% identity), but less with the mammalian AdoMetDC (30–33% identity) [10].
Till now, various compounds have been investigated to inhibit AdoMetDC in in vitro studies such as carbonimidic dihydrazide (CGP40215) and 5-(((Z)-4-amino-2-butenyl)methylamino)- 5-deoxyadenosine (MDL73811) [13,18-20]. Several other trypanocidal drugs, including methylglyoxyl bis-guanylhydrazone, berenil, and pentamidine, were also investigated as AdoMetDC inhibitors [8,9]. Pentamidine, the second most frequently prescribed drug, is also used for the treatment of visceral leishmaniasis [4,21]. It has been suggested that pentamidine and berenil may have other targets within the cell, and the observed toxicity to the host may also be due to their potential lack of specificity toward a single target [22,23].
The antileishmanial activity of different plant extracts has been observed in various reports, although the precise mechanisms through which these extracts combat the disease remain challenging to understand due to the complex composition of phytochemicals present in crude plant extracts. Phytochemicals, or natural products, have always played a significant role in the treatment of various diseases [24,25], exhibiting a wide range of pharmacological properties such as antimicrobial, antioxidant, anticarcinogenic, and anti-inflammatory activities [26,27]. In this study, our objective is to identify specific inhibitors of LdAdoMetDC from medicinal plants that have been previously reported for their antileishmanial activities. By focusing on these specific inhibitors, we aim to shed light on the molecular mechanisms underlying the observed antileishmanial effects of these medicinal plants.
2. MATERIALS AND METHODS
2.1. Homology Modeling and Ligand Binding Site Identification
Since the crystal structure of LdAdoMetDC is not available, its protein sequence (Accession no. TPP45862.1, Length: 382 amino acids) has been retrieved from the NCBI protein database. The suitable templates for homology modeling have been searched by blasting the LdAdoMetDC protein sequence against the protein data bank using PSI-BLAST [28]. Crystal structure of Trypanosoma brucei AdoMetDC (TbAdoMetDC) at 2.42 Å resolutions (5TVF_A and 5TVF-B) has been selected as templates based on the query coverage (21% and 74%, respectively) and identity of 57.14% and 65.14%, respectively, with the sequence of LdAdoMetDC. Together, the templates (5TVF_A and 5TVF_B) provide the 95% query coverage [Figure 1]. Twenty-five homology models of LdAdoMetDC have been prepared using Modeller 10.1, through a multi-template homology modeling approach. Model that exhibited the lowest discrete optimized protein energy (DOPE score: −42092.67969) has been taken for further optimization [29]. Energy minimization of this model has been done using YASARA web server which performs energy minimization of protein models in explicit solvent using its own developed optimized force field [30]. The predicted model has been further evaluated for its quality using PROCHECK available at Structure Analysis and Verification server (SAVES v6.0) and ProSA-web server [31-34]. In addition, root mean square deviation (RMSD) has been obtained by structurally aligning the predicted model with the template using PyMOL. The ligand binding site of LdAdoMetDC has been predicted by superimposing the modeled structure of LdAdoMetDC with ligand-bound (CGP40215) crystal structure of TbAdoMetDC.
 | Figure 1: Multiple sequence alignment of LdAdoMetDC with template (Trypanosoma brucei AdoMetDC) 5TVF_Chain-A and 5TVF_Chain-B. [Click here to view] |
2.2. Selection of Plants and Phytochemicals
Based on the available literature of folk medicine, 155 medicinal plants with antileishmanial activities have been selected. One thousand eight hundred and sixty-seven phytochemicals of these plants have been downloaded in PDB file format from IMPPAT database which is one of the most extensive resources available till date [35]. Known inhibitors of LdAdoMetDC reported in different literatures have been obtained in PDB file format using UCSF Chimera v1.15 after being retrieved in SDF file format from the PubChem [9,18,19,36-38]. A library has been prepared collecting PDB files of phytochemicals and known inhibitors of AdoMetDC.
2.3. Molecular Docking
Molecular docking of all the phytochemicals along with the known inhibitors has been performed using a graphical user interface (GUI) based tool Raccoon 1.0. It utilizes AutoDock Tools (ADT) 1.5.6 for preparation of parameter files and AutoDock 4.2.1 as a tool for docking [39,40]. Raccoon automatically processed ligand libraries and generated PDBQT input files after adding polar hydrogens and assigning Gasteiger charges to all the small molecules. Using ADT, the target macromolecule has been prepared separately by saving in PDBQT file format after the allocation of Gasteiger charges. During the docking, rotatable bonds of all the phytochemicals have been considered as rotatable and the target macromolecule has been considered as rigid. The configuration files for the grid parameters and docking parameters have been generated using the ADT. Raccoon automatically generates grid maps for each of the ligands. Grid box size of 100 × 100 × 100 Å with 0.375 Å spacing has been selected that covers all the identified binding site residues. Further, molecular docking has been performed using Lamarckian genetic algorithm and empirical-free energy functions. The process began with an initial 150 randomly placed individual’s population, followed by a maximum energy evaluation of 2,500,000. The crossover rate and mutation rate has been set at 0.80 and 0.02, respectively. For each phytochemicals, 10 separate docking runs have been carried out remaining all the value of parameters as default. Based on binding free energy (ΔG), top 15 phytochemicals have been selected and carried forward along with eight known inhibitors of LdAdoMetDC for 100 separate docking runs remaining all the parameters same as mentioned above. The LdAdoMetDC-phytochemical complexes and LdAdoMetDC-known inhibitor complexes with the lowest DG value from the largest cluster have been written in PDBQT format and converted to PDB file format using PyMOL. Further, these complexes have been analyzed using PyMOL for possible polar and hydrophobic interactions. All the docking studies have been performed at Intel (R) Core (TM) i7-3770 CPU (3.40 GHz) with Linux-based operating system Ubuntu 18.04 LTS.
2.4. MD Simulations
MD simulations have been performed to explore the binding stability and dynamic behavior of unbound LdAdoMetDC, LdAdoMetDC-phytochemicals complexes, and LdAdoMetDC known inhibitor complex using a freely available online server WebGRO [41,42]. The GROMACS software is utilized by this server to perform molecular dynamics (MD) simulations of protein and protein-ligand complexes [43,44]. The topologies and parameters for the macromolecule have been generated using GROMOS96 54a7 force field and for ligands using PRODRG server [45,46]. Each unbound macromolecule/docked complex has been put inside a cubic simulation box with edges spaced apart by a factor of 1 Å distance. During the solvation process, simple point charge water model has been selected and the system has been neutralized by adding counter ions (NaCl) to bring a molarity of 0.15 M. Further, the energy of each unbound macromolecule/docked complex has been minimized using steepest descent integrator at every 50,000 steps followed by 50,000 steps of equilibration both in NVT (constant number of particles, temperature, and volume) and NPT (constant number of particles, pressure, and temperature) ensemble [47]. The temperature and pressure of each system has been controlled by a Berendsen thermostat and a Parrinello-Rahman barostat, which have been set at 300 K and 1 bar, respectively [44,48]. Further, using Leap-frog integrator, each system has been simulated for 100 ns, and the number of frames generated per simulation were 5,000 [49]. The GROMACS analytic methods have been employed to determine the RMSD, root-mean-square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and hydrogen bonds from the generated trajectories [43]. Visual Molecular Dynamics (VMD) software has been used to visualize the trajectories and Grace v5.1.25 has been used to create the graphs [50,51].
2.5. Binding Free Energy Calculation
The molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method (g_mmpbsa) has been used for calculating the binding free energy by reading the MD trajectories. It calculates binding energy components as well as residue wise energy contributions. MM-PBSA method has gained recognition for its ability to more precisely predict the free energies of ligand binding, in contrast to other computational methods such as docking [52,53].
3. RESULTS AND DISCUSSION
3.1. Homology Model of LdAdoMetDC and Its Ligand Binding Site Residues
The predicted model [Figure 2a] has been evaluated for its quality using various programs available at SAVES v6.0 server. The Ramachandran plot analysis revealed that the largest proportion of residues (89.7%) were located in the most preferred regions. In addition, 8.7% of the residues were located in other allowed regions, while 1.0% and 0.6% were found in generously allowed regions and disallowed regions, respectively [Figure 2b]. Furthermore, an interactive web service, ProSA-web, has been used for detecting errors in three-dimensional protein structures. The model’s calculated z-score was −7.45, as shown in [Figure 2c], and it has been predicted to be of X-ray crystallographic structure quality. An additional approach to assess the quality of the LdAdoMetDC model locally involved generating a plot depicting the knowledge-based energies in relation to the position of the amino acid residues. The majority of the residues had negative energy values predicted [Figure 2d]. The RMSD of the template model was 0.066 Å. According to the analysis, the modeled structure was of high quality. The binding site residues of LdAdoMetDC were Phe-30, Phe-32, Leu-87, Thr-88, Glu-89, Cys-104, Phe-248, Pro-250, Cys-251, Gly-252, Tyr-253, Ser-254, His-267, Ile-268, Thr-269, Pro-270, and Glu-271 [Figure 3].
 | Figure 2: Structure validation of predicted model using SAVESv6.0 and ProSA-web server. (a) Predicted LdAdoMetDC model, (b) Ramachandran plot of LdAdoMetDC obtained through PROCHECK representing the percentages of residues in the most favored regions (89.7%), additionally allowed regions (8.7%), generously allowed regions (1.0%), and disallowed regions (0.6%) (c) Z-score plot obtained throuth ProSA-web, depicting overall quality of the model (d) Knowledge-based energy plot obtained through ProSA-web depicting local model quality. [Click here to view] |
 | Figure 3: Prediction of active site amino acid residues (yellow) by superimposing LdAdoMetDC model with the template around CID9576798 (CGP40215), inhibitor of TbAdoMetDC (5TVF). [Click here to view] |
3.2. Molecular Docking Studies
Binding free energy (ΔG), possible polar and hydrophobic interactions of top 15 docked phytochemicals, and 8 known inhibitors docked to LdAdoMetDC are shown in Tables 1 and 2, respectively. The results demonstrated that CID5488537 (Fagopyrine), CID442630 (Carpaine), and CID44558930 (Anabsinthin) were the molecules with the lowest ΔG of −9.42, −8.96 and −8.66 kcal/mol, respectively, among all the docked phytochemicals. However, among the known inhibitors CID9576798 (CGP-40215) docked with the lowest ΔG of −6.75 kcal/mol. It has been observed that the binding modes of these phytochemicals were similar to the known inhibitor [Figure 4a]. The binding mode of CID5488537 showed that the hydroxyl group at one aromatic ring made polar interaction with the side chain of Glu-89. However, on the other side of the molecule, two hydroxyl groups and one carbonyl group present separately at three aromatic rings were in polar interaction range with the side chain of Glu-271. Furthermore, the nitrogen atom of the piperidine group was in polar interaction with the side chain of Glu-249 [Figure 4b]. The binding mode of CID442630 revealed that the nitrogen at one of the azatricyclo ring of the molecule established polar interaction with the side chain of Glu-89 [Figure 4c]. When CID44558930 docked to LdAdoMetDC, the hydroxyl group at the cycloheptane ring made polar interaction with the side chain of Glu-249 [Figure 4d]. Among the docked known inhibitors, CID9576798 binding mode studies revealed that the -NH2 group on one side of the molecule established polar interaction with the side chain of Gly-29, while the nitrogen atom on the other side of the molecule made polar interaction with Leu-87. The -NH2 group and a nitrogen atom in the middle of the molecule was in polar interaction range with Glu-249 [Figure 4e]. The majority of docked phytochemicals formed polar interactions with amino acid residues such as Glu-89, Glu-249, His-267, and Ser-254, whereas hydrophobic residues at the binding site included Cys-104, Cys-251, Gly-252, Ile-268, Phe-30, Phe-32, Phe-248, and Pro-250.
Table 1: Binding free energy (ΔG) estimated with AutoDock 4.2 and interaction of phytochemicals with LdAdoMetDC predicted by PyMOL.
| Phytochemicals | ΔG (kcal/mol) | Putative Polar Interactions | Hydrophobic residues in 4Å region |
|---|---|---|---|
| CID5488537 (Fagopyrine) | −9.42 | Glu-89, Glu-271, Glu-249 | Gly-29, Phe-30, Phe-32, Phe-217, Phe-248, Pro-250 |
| CID442630 (Carpaine) | −8.96 | Glu-89 | Gly-252, Ile-268, Phe-30, Phe-32, Phe-248, Pro-250 |
| CID44558930 (Anabsinthin) | −8.66 | Glu-249 | Gly-252, Ile-268, Phe-32, Phe-217, Phe-248 |
| CID101316729 ((1R,3aR,5aR,5bS,7aS,9S,11aS,11bR,13aS,13bR)-3a, 5a, 7a, 11b, 13a-pentamethyl-8-methylidene-1- propan-2-yl-2,3,4,5,5b, 6,7,9,10,11,11a, 12,13,13b-tetradecahydro- 1H-cyclopenta[a] chrysen-9-ol) | −8.42 | His-267, Ser-254 | Gly-252, Ile-268, Phe-30, Phe-32, Phe-248, Pro-250 |
| CID14034468 (Isomangiferolic acid) | −8.28 | Gly-29, His-267, Leu-89 | Gly-252, Leu-87, Phe-30, Phe-32, Phe-248, Pro-250 |
| CID3035446 (Sarsaponin) | −8.28 | Ile-268, Ser-254 | Gly-252, Phe-30, Phe-32, Phe-248 |
| CID14034474 (Mangiferonic acid) | −8.27 | Ile-102, Leu-87, Ser-91, Val-86 | Gly-252, Phe-30, Phe-32, Phe-248, Pro-250 |
| CID21679023 (Withanolide G) | −7.96 | Glu-89, Glu-249, Leu-87 | Phe-30, Phe-32, Phe-217, Phe-248, Pro-250 |
| CID23266155 (27-Deoxywithaferin A) | −7.99 | Gly-29, Ser- 254 | Gly-252, Ile-268, Phe-30, Phe-32, Phe-248, Pro-250 |
| CID10895555 (Dammarenediol II) | −8.24 | Glu-89, Glu-271 | Gly-252, Ile-268, Phe-32, Phe-217, Phe-248, Pro-250, Pro-270 |
| CID160482 (Lophenol) | −8.22 | His-267, Ser-254, Tyr-253 | Gly-29, Gly-252, Ile-268, Phe-30, Phe-32, Phe-248 |
| CID5283652 (24-methylcholesta-5,23E-dien-3beta-ol) | −8.03 | His-267, Ser- 254 | Gly-29, Gly-252, Ile-268, Phe-30, Phe-32, Phe-248, Pro-250, |
| CID101286241 ((2R,6R)-6-[(1S,3R,6S,8R,11S,12S,15R,16R) - 6-hydroxy-7,7,12,16-tetramethyl-15- pentacyclo[9.7.0.01,3.03,8.012,16]octadecanyl]-2-methyl-3- methylideneheptanoic acid) | −8.04 | Leu-87, His-267 | Gly-252, Ile-268, Phe-32, Phe-248, Phe_30, Pro-250 |
| CID21159864 (A-spinasterone) | −8.00 | His-267 | Gly-29, Gly-252, Ile-268, Phe-30, Phe-32, Phe-248, Pro-250 |
| CID10478550 (Arnidenediol) | −7.98 | Glu-249, Gly-29 | Phe-30, Phe-32, Phe-217, Phe-248, Pro-250 |
Table 2: Free energy of binding (ΔG) estimated with AutoDock 4.2 and interaction of known inhibitors with LdAdoMetDC predicted by PyMOL.
| Known inhibitors | ΔG (kcal/mol) | Putative polar interactions | Hydrophobic residues in 4Å region |
|---|---|---|---|
| CID9576798 (CGP-40215; Carbonimidic dihydrazide) | −6.75 | Glu-249, Gly-29, Leu-87 | Gly-252, Ile-268, Phe-30, Phe-32, Phe-217, Phe-248, Pro-250 |
| CID9576789 (CGP-48664; Sardomozide) | −6.34 | Glu-271, Leu-87 | Gly-252, Ile-268, Phe-32, Phe-248 |
| CID6436013 (MDL-73811) | −5.84 | Glu-89, Glu-249, Glu-271, Leu-87 | Gly-252, His-267, Ile-268, Phe-32, Phe-248, Pro-250 |
| CID2354 (Berenil) | −5.95 | Glu-271, Ser-254, Thr-107 | Ala-71, Cys-104, Cys-251, Gly-252, Ile-268, Leu-87, Phe-32, Phe-248, Pro-270 |
| CID5351154 (Mitoguazone) | −5.31 | Cys-251, Glu-249, Glu-271, Ile-268 | Gly-252, Phe-32, Phe-248, Pro-250, Pro-270 |
| CID4735 (Pentamidine) | −4.26 | Glu-89, Glu-249, Glu-271, Ile-268 | Gly-252, Phe-32, Phe-217, Phe-248, Pro-250, Ser-254 |
| CID122092 (MHZPA) | −4.54 | Cys-251, Glu-89, Glu-249, Glu-271 | Gly-252, Ile-268, Phe-32, Phe-248, Pro-250, Pro-27 |
| CID65482 (Sinefungin) | −3.56 | Cys-251, Glu-89, Glu-249, Ile-268, Thr-269 | Gly-252, Phe-32, Phe-217, Phe-248, Pro-250, Pro-270 |
 | Figure 4: Binding mode analysis (a) Superimposed binding modes of the top three phytochemicals, CID5488537 (yellow), CID442630 (cyan), and CID44558930 (blue) and one known inhibitor, CID9576798 (orange) at the active site of LdAdoMetDC (b) CID5488537 showing polar contacts with Glu-89, Glu-271 and Glu-249 (c) CID442630 showing polar contacts with Glu-89 (d) CID44558930 showing polar contacts with Glu-249 (e) CID9576798 (orange) showing polar contacts with Cys-251, Glu-249, Gly-29, and Leu-87. [Click here to view] |
3.3. MD Simulations of Proteins and Protein–Ligand Complexes
100 ns MD simulations have been used to investigate the dynamic behavior of LdAdoMetDC and LdAdoMetDC-ligand complexes. RMSD, RMSF, Rg, hydrogen bond, and SASA have been calculated for the LdAdoMetDC-ligand complexes.
3.3.1. RMSD analysis
RMSD is a commonly used metric to assess changes between the reference structure and the structures sampled throughout the simulation. MD simulations have been done to make sure the docked LdAdoMetDC-phytochemical complexes which were dynamically stable or not. The MD simulation has also been run for the known inhibitor complexes for comparative analysis. The RMSD time trajectory reflects the variation between a protein and ligand structure with a reference structure over time [57]. The backbone of the docked protein and ligands structures have been used as a reference to generate RMSD graphs with respect to 100 ns production run time. Figure 5 shows the graphically superimposed time-dependent RMSD of LdAdoMetDC alone and LdAdoMetDC-ligand complexes. The RMSD values of all four structures increased gradually from 0 to 40 ns with minor fluctuations. After 40 ns LdAdoMetDC-CID9576798, LdAdoMetDC-CID442630 and LdAdoMetDC-CID44558930 complexes converge and attained stability with an average RMSD of 0.2823, 0.2907, and 0.2865, respectively, at the end of the 100 ns simulation run with no significant fluctuations, while LdAdoMetDC- CID5488537 complex has been observed to rise steadily after 40 ns and attained stability with an average RMSD of 0.3915 at the end of 100 ns MD run. However, LdAdoMetDC alone stabilized after 23 ns of simulation time and remained in this state for the rest of the simulation with average RMSD value of 0.2980 nm. No significant variations have been observed in the relative RMSD values of various complexes except the LdAdoMetDC-CID5488537 complex. The differences in the RMSD value of LdAdoMetDC-CID5488537 complex relative to the RMSD of other three complexes demonstrated that the ligand binding had an impact on the corresponding LdAdoMetDC structure.
 | Figure 5: Root mean square deviation plot of LdAdoMetDC alone (black), LdAdoMetDC-CID5488537 (green), LdAdoMetDC-CID442630 (blue), LdAdoMetDC-CID44558930 (orange), and LdAdoMetDC-CID9576798 (red) complexes generated over a 100 ns molecular dynamics simulation (*known inhibitor complex). [Click here to view] |
3.3.2. RMSF analysis
RMSF indicates the pliability of various segments of a protein and correlates it with the B-factors observed in crystallography. The RMSF trajectories of docked complexes have been studied. The assessment of the stability profile has been utilized to examine the amino acid residues that play a role in the intricate changes in structure. Greater variations are indicated by higher RMSF values. The superimposed RMSF value per residue for docked complexes is shown in Figure 6. The average RMSF values for LdAdoMetDC alone, LdAdoMetDC-CID5488537, LdAdoMetDC-CID442630, LdAdoMetDC-CID44558930, and LdAdoMetDC-CID9576798 complex were 0.1104 nm, 0.1710 nm, 0.1610 nm, 0.1439 nm, and 0.1560 nm, respectively. The RMSF values of the active site residues for the LdAdoMetDC-phytochemical complexes showed similar fluctuating patterns with those of the known inhibitor.
 | Figure 6: Root mean square fluctions trajectories of residues for LdAdoMetDC alone (black), LdAdoMetDC-CID5488537 (green), LdAdoMetDC-CID442630 (blue), LdAdoMetDC-CID44558930 (orange), and LdAdoMetDC-CID9576798 (red) complexes (*known inhibitor complex). [Click here to view] |
3.3.3. Rg analysis
Using the Rg value, the compactness and structural alterations of the docked complexes have been evaluated. It is a measure of determining the mass of atoms in relation to the center of mass in a complex protein. Proteins that are folded show tight packing, whereas proteins that are unfolded show loose packing, less stable conformation, and larger values for the Rg. It has been observed that throughout the MD simulation, the Rg value of LdAdoMetDC alone and LdAdoMetDC-ligand complexes remained mostly stable [Figure 7]. This suggests that the binding of phytochemicals and known inhibitor do not cause any significant structural changes to LdAdoMetDC and that it has remained structurally stable in its complex state with these ligands. LdAdoMetDC alone, LdAdoMetDC-CID9576798, LdAdoMetDC-CID5488537, LdAdoMetDC-CID442630, and LdAdoMetDC-CID44558930 complexes have shown an average Rg value of 1.8902, 1.9664 nm, 1.9621 nm, 1.9521 nm, and 1.9416 nm, respectively, for 100 ns of simulation duration, which are quite close and do not exhibit any significant differences.
 | Figure 7: Radius of gyration (Rg) trajectories versus time graph for LdAdoMetDC alone (black), LdAdoMetDC-CID5488537 (green), LdAdoMetDC-CID442630 (blue), LdAdoMetDC-CID44558930 (orange), and LdAdoMetDC-CID9576798 (red) complexes (*known inhibitor complex). [Click here to view] |
3.3.4. Hydrogen bond analysis
In a subsequent analysis, the intermolecular hydrogen bond formation has been evaluated for each of the four (one known and three phytochemicals) complexes during the 100 ns simulation run have been depicted in Figure 8. The intermolecular hydrogen bonds that were seen to form in the LdAdoMetDC-CID9576798 complex, on average were four, with a maximum of six. Moreover, the intermolecular hydrogen bonds in LdAdoMetDC-CID5488537, LdAdoMetDC-CID442630 and LdAdoMetDC-CID44558930 complexes showed an average of two, one and two, while a maximum of six, two, and four, respectively.
 | Figure 8: Hydrogen bonds trajectories versus time graph for LdAdoMetDC-CID5488537 (green), LdAdoMetDC-CID442630 (blue), LdAdoMetDC-CID44558930 (orange), and LdAdoMetDC-CID9576798 (red) complexes (*known inhibitor complex). [Click here to view] |
3.3.5. SASA analysis
The change in the SASA value for the LdAdoMetDC complexes during the 100 ns simulation time is shown in Figure 9. The superimposed plot showed that the value of SASA for the phytochemical complexes, LdAdoMetDC-CID5488537 (162.122 ± 7.668 nm2), LdAdoMetDC-CID442630 (158.235 ± 5.568 nm2), and LdAdoMetDC-CID44558930 (158.655 ± 5.988 nm2) was lower than the value for the known inhibitor complex, LdAdoMetDC-CID9576798 (165.378 ± 5.042 nm2) at the end of 100 ns of MD run. All the complexes showed minor fluctuations throughout the 100 ns MD run. Binding of phytochemicals and known inhibitor to LdAdoMetDC caused insignificant changes in the protein during the course of the production period of MD simulations. The results of SASA indicated that phytochemical complexes, LdAdoMetDC-CID442630 followed by LdAdoMetDC-CID44558930 and LdAdoMetDC-CID5488537, were consistently more stable than known inhibitor complex, LdAdoMetDC-CID9576798.
 | Figure 9: Superimposed SASA trajectories versus time graph for LdAdoMetDC alone (black), LdAdoMetDC-CID5488537 (green), LdAdoMetDC-CID442630 (blue), LdAdoMetDC-CID44558930 (orange), and LdAdoMetDC-CID9576798 (red) complexes (*known inhibitor complex). [Click here to view] |
3.3.6. Binding free energy calculation
The MD trajectories have been utilized to determine the binding free energy of the simulated complexes, to reaffirm the inhibitor’s affinity that has been anticipated by the docking studies. The MM-PBSA method has been employed to calculate the overall non-polar, polar, and non-bonded interaction energies (which includes van der Waals and electrostatic interactions) of each complex [Table 3]. The computed binding free energies of the phytochemicals complexes, LdAdoMetDC-CID5488537, LdAdoMetDC-CID442630 and LdAdoMetDC-CID44558930 were -233.853, -131.805 and -149.506 kJ/mol respectively. However, the binding free energy of the LdAdoMetDC-CID9576798 complex was −125.695 kJ/mol. The binding energies determined from MD simulations are consistent with the docking results. It has been observed that the total binding free energies of the complexes were greatly impacted by Van Der Waals and electrostatic interactions.
Table 3: Components of binding free Energy (kJ/mol) for the top three phytochemical complexes and one known inhibitor complex.
| Energies in (kJ/mol) | LdAdoMetDC-CID5488537 | LdAdoMetDC-CID442630 | LdAdoMetDC-CID44558930 | LdAdoMetDC-CID9576798* |
|---|---|---|---|---|
| Van Der Waal energy | −257.548 | −155.105 | −179.090 | −148.031 |
| Electrostattic energy | −7.737 | −25.784 | −26.886 | −3.168 |
| Polar solvation energy | 54.104 | 69.327 | 73.699 | 41.412 |
| Nonpolar solvation energy | −22.669 | −20.251 | −17.218 | −15.969 |
| Binding free energy total (ΔG) | −233.853 | −131.805 | −149.506 | −125.695 |
* Known inhibitor
3.3.7. Analysis of residue-wise binding energy contribution
Contributions of binding energy per residue have been estimated for all the complexes using MM-PBSA method [58]. Residues that play a crucial role in the binding of a ligand to a protein, by contributing binding free energy of approximately ±5 kJ/mol, can be considered as key residues [59]. The binding free energy contributions by active site residues in the complexes are shown in Figure 10. For the LdAdoMetDC-CID9576798 complex, it has been found that Phe-32, Glu-89, and Phe-248 contributed energies that were above the threshold of ±5 kJ/mol, with the value of −6.6949, −5.2507, and −6.0216 kJ/mol, respectively. In the LdAdoMetDC-CID5488537 complex, Phe-32 and Phe-248 contributed energies of −14.6285 and −18.1273 kJ/mol, respectively. For the LdAdoMetDC-CID442630 complex, Phe32, Glu-89, and Glu-271 were shown to contribute energy above the threshold of ±5 kJ/mol with values of −7.0677, −12.591, and −5.1948 kJ/mol, respectively. In the LdAdoMetDC-CID44558930 complex, Glu-89 and Phe-248 were the residues that contribute the energy over the threshold of ±5 kJ/mol, with an energy value of −11.1276 and −8.6428 kJ/mol.
 | Figure 10: Binding free energy (ΔG) contribution from the active site residues for LdAdoMetDC -CID5488537 (green), LdAdoMetDC-CID442630 (blue), LdAdoMetDC-CID44558930 (orange), and LdAdoMetDC-CID9576798 (red) complexes (*known inhibitor complex). [Click here to view] |
3.4. Drug-Likeness and ADME Analysis
Physicochemical properties have been summed up using SwissADME in Table 4 to assess the drug-likeness of the studied phytochemicals. CID442630 and CID44558930 followed all the Lipinski’s Rule parameters as well as Veber rule parameters whilst CID5488537 violated two parameters of Lipinski rule as well as one parameter of Veber rule as it has molecular weight of 670.71 g/mol, 8 hydrogen bond donors, and TPSA of 179.58 Ų. Thus, CID5488537 has been predicted to have poor bioavailability as well as cell membrane permeability.
Table 4: Physicochemical properties of top docked phytochemicals predicted by SwissADME.
| Phytochemicals | LogP | Molecular weight g/mol | H-bond acceptors | H-bond donors | Rotatable bonds | TPSA (Ų) |
|---|---|---|---|---|---|---|
| CID5488537 | 1.61 | 670.71 | 10 | 8 | 2 | 179.58 |
| CID442630 | 3.75 | 478.71 | 6 | 2 | 0 | 76.66 |
| CID44558930 | 4.04 | 496.64 | 6 | 1 | 0 | 82.06 |
4. CONCLUSIONS
In the presented work, the interactions of phytochemicals from the plants with antileishmanial activities have been explored by molecular docking and MD simulation studies. It has been found that CID5488537, CID442630, and CID44558930 exhibit the best binding affinities in molecular docking studies. CID5488537, CID442630, and CID44558930 belong to medicinal plants, Fagopyrum esculentum (Buckwheat), Carica papaya (Papaya), and Artemisia absinthium (Wormwood), respectively [60-62]. MD simulations have been performed to analyze the stability of LdAdoMetDC-phytochemical complexes and compared with that of LdAdoMetDC known inhibitor complex. MD results demonstrated stable RMSD, RMSF, Rg, and SASA for docked LdAdoMetDC complexes. Further, hydrogen bond analysis has been employed to evaluate the protein-ligand interactions, and binding free energy has been calculated to determine the binding affinity based on MD trajectories. CID5488537, CID442630, and CID44558930 exhibited better binding free energy in comparison to known inhibitor (CID9576798) calculated by MM-PBSA method. It has also been observed that Van Der Waals interactions made significant contributions in binding free energies of the docked complexes. Furthermore, these molecules might function as potential lead molecules for the development of potential LdAdoMetDC inhibitors.
5. AUTHORS’ CONTRIBUTIONS
PKA: Data collection, compilation, drafting, and art work of the manuscript. KB: Manuscript Editing. AKS: Manuscript moderation. AK: Primary investigator of the work presented.
6. FUNDING
There is no funding to report.
7. CONFLICTS OF INTEREST
The authors report no financial or any other conflicts of interest in this work.
8. ETHICAL APPROVALS
This study does not involve experiments on animals or human subjects.
9. DATA AVAILABILITY
All the data is available with the authors and shall be provided upon request.
10. PUBLISHER’S NOTE
This journal remains neutral with regard to jurisdictional claims in published institutional affiliation.
REFERENCES
1. Scarpini S, Dondi A, Totaro C, Biagi C, Melchionda F, Zama D, et al. Visceral leishmaniasis:Epidemiology, diagnosis, and treatment regimens in different geographical areas with a focus on pediatrics. Microorganisms 2022;10:1887. [CrossRef]
2. Sunter J, Gull K. Shape, form, function and Leishmania pathogenicity:From textbook descriptions to biological understanding. Open Biol 2017;7:170165. [CrossRef]
3. Moreira PO, Nogueira PM, Monte-Neto RL. Next-generation leishmanization:Revisiting molecular targets for selecting genetically engineered live-attenuated Leishmania. Microorganisms 2023;11:1043. [CrossRef]
4. Verdan M, Taveira I, Lima F, Abreu F, Nico D. Drugs and nanoformulations for the management of Leishmania infection:A patent and literature review (2015-2022). Expert Opin Ther Pat 2023;33:137-50. [CrossRef]
5. Pawar S, Kumawat MK, Kundu M, Kumar K. Synthetic and medicinal perspective of antileishmanial agents:An overview. J Mol Struct 2023;1271:133977. [CrossRef]
6. Majumder N, Banerjee A, Saha S. A review on new natural and synthetic anti-leishmanial chemotherapeutic agents and current perspective of treatment approaches. Acta Trop 2023;240:106846. [CrossRef]
7. Tali MB, Kamdem BP, Tchouankeu JC, Boyom FF. Current developments on the antimalarial, antileishmanial, and antitrypanosomal potential and mechanisms of action of Terminalia spp. S Afr J Bot 2023;156:309-33. [CrossRef]
8. Carter NS, Kawasaki Y, Nahata SS, Elikaee S, Rajab S, Salam L, et al. Polyamine metabolism in Leishmania parasites:A promising therapeutic target. Med Sci (Basel) 2022;10:24. [CrossRef]
9. Roberts SC, Scott J, Gasteier JE, Jiang Y, Brooks B, Jardim A, et al. S-adenosylmethionine decarboxylase from Leishmania donovani. Molecular, genetic, and biochemical characterization of null mutants and overproducers. J Biol Chem 2002;277:5902-9. [CrossRef]
10. Persson L. Ornithine decarboxylase and S-adenosylmethionine decarboxylase in trypanosomatids. Biochem Soc Trans 2007;35:314-7. [CrossRef]
11. Soni M, Pratap JV. Development of novel anti-leishmanials:The case for structure-based approaches. Pathogens 2022;11:950. [CrossRef]
12. Colotti G, Ilari A. Polyamine metabolism in Leishmania:From arginine to trypanothione. Amino Acids 2010;40:269-85. [CrossRef]
13. Singh SP, Agnihotri P, Pratap JV. Characterization of a novel putative S-adenosylmethionine decarboxylase-like protein from Leishmania donovani. PLoS One 2013;8:e65912. [CrossRef]
14. Gupta D, Singh PK, Yadav PK, Narender T, Patil UK, Jain SK, et al. Emerging strategies and challenges of molecular therapeutics in antileishmanial drug development. Int Immunopharmacol 2023;115:109649. [CrossRef]
15. Gilroy C, Olenyik T, Roberts SC, Ullman B. Spermidine synthase is required for virulence of Leishmania donovani. Infect Immun 2011;79:2764-9. [CrossRef]
16. Taladriz S, Ramiro MJ, Hanke T, Larraga V. S-adenosylmethionine decarboxylase from Leishmania infantum promastigotes:Molecular cloning and differential expression. Parasitol Res 2002;88:421-6. [CrossRef]
17. Mishra AK, Agnihotri P, Srivastava VK, Pratap JV. Novel protein-protein interaction between spermidine synthase and S-adenosylmethionine decarboxylase from Leishmania donovani. Biochem Biophys Res Commun 2015;456:637-42. [CrossRef]
18. Mukhopadhyay R, Kapoor P, Madhubala R. Antileishmanial effect of a potent S-adenosylmethionine decarboxylase inhibitor:CGP 40215A. Pharmacol Res 1996;33:67-70. [CrossRef]
19. Mukhopadhyay R, Madhubala R. Antileishmanial activity of berenil and methylglyoxal bis (guanylhydrazone) and its correlation with S-adenosylmethionine decarboxylase and polyamines. Int J Biochem Cell Biol 1995;27:55-9. [CrossRef]
20. Pardali V, Giannakopoulou E, Balourdas DI, Myrianthopoulos V, Taylor MC, Šekutor M, et al. Lipophilic guanylhydrazone analogues as promising trypanocidal agents:An extended SAR study. Curr Pharm Des 2020;26:838-66. [CrossRef]
21. Akbari M, Oryan A, Hatam G. Application of nanotechnology in treatment of leishmaniasis:A review. Acta Trop 2017;172:86-90. [CrossRef]
22. Basselin M, Robert-Gero M. Alterations in membrane fluidity, lipid metabolism, mitochondrial activity, and lipophosphoglycan expression in pentamidine-resistant Leishmania. Parasitol Res 1997;84:78-83. [CrossRef]
23. Singh R, Siddiqui KA, Valenzuela MS, Majumder HK. Kinetoplast DNA minicircle binding proteins in a Leishmania spp:Interference of protein DNA interaction by berenil. Indian J Biochem Biophys 1995;32:437-41.
24. Islam MT, Sarkar C, El-Kersh DM, Jamaddar S, Uddin SJ, Shilpi JA, et al. Natural products and their derivatives against Coronavirus:A review of the non-clinical and pre-clinical data. Phytother Res 2020;34:2471-92. [CrossRef]
25. Yuan H, Ma Q, Ye L, Piao G. The traditional medicine and modern medicine from natural products. Molecules 2016;21:559. [CrossRef]
26. Bourais I, Elmarrkechy S, Taha D, Mourabit Y, Bouyahya A, El Yadini M, et al. A review on medicinal uses, nutritional value, and antimicrobial, antioxidant, anti-inflammatory, antidiabetic, and anticancer potential related to bioactive compounds of J. regia. Food Rev Int 2022;2022:1-51. [CrossRef]
27. Ahmad B, Rehman MU, Amin I, Arif A, Rasool S, Bhat SA, et al. A review on pharmacological properties of zingerone (4-(4-hydroxy-3-methoxyphenyl)-2-butanone). ScientificWorldJournal 2015;2015:816364. [CrossRef]
28. Margelevicius M, Venclovas C. PSI-BLAST-ISS:An intermediate sequence search tool for estimation of the position-specific alignment reliability. BMC Bioinformatics 2005;6:185. [CrossRef]
29. Webb B, Sali A. Comparative protein structure modeling using MODELLER. Curr Protoc Bioinformatics 2016;54:5.6.1-37. [CrossRef]
30. Krieger E, Koraimann G, Vriend G. Increasing the precision of comparative models with YASARA NOVA--a self-parameterizing force field. Proteins 2002;47:393-402. [CrossRef]
31. Laskowski RA, MacArthur MW, Moss DS, Thornton JM. PROCHECK:A program to check the stereochemical quality of protein structures. J Appl Crystallogr 1993;26:283-91. [CrossRef]
32. Colovos C, Yeates TO. Verification of protein structures:Patterns of nonbonded atomic interactions. Protein Sci 1993;2:1511-9. [CrossRef]
33. Wiederstein M, Sippl MJ. ProSA-web:Interactive web service for the recognition of errors in three-dimensional structures of proteins. Nucleic Acids Res 2007;35:W407-10. [CrossRef]
34. Sippl MJ. Recognition of errors in three-dimensional structures of proteins. Proteins 1993;17:355-62. [CrossRef]
35. Mohanraj K, Karthikeyan BS, Vivek-Ananth RP, Chand RP, Aparna SR, Mangalapandi P, et al. IMPPAT:A curated database of Indian medicinal plants, phytochemistry and therapeutics. Sci Rep 2018;8:4329. [CrossRef]
36. Millward MJ, Joshua A, Kefford R, Aamdal S, Thomson D, Hersey P, et al. Multi-centre phase II trial of the polyamine synthesis inhibitor SAM486A (CGP4?) in patients with metastatic melanoma. Invest New Drugs 2005;23:253-6. [CrossRef]
37. Phelouzat MA, Lawrence F, Moulay L, Borot C, Schaeverbeke J, Schaeverbeke M, et al. Leishmania donovani:Antagonistic effect of S-adenosyl methionine on ultrastructural changes and growth inhibition induced by sinefungin. Exp Parasitol 1992;74:177-87. [CrossRef]
38. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera--a visualization system for exploratory research and analysis. J Comput Chem 2004;25:1605-12. [CrossRef]
39. Forli S, Huey R, Pique ME, Sanner MF, Goodsell DS, Olson AJ. Computational protein-ligand docking and virtual drug screening with the AutoDock suite. Nat Protoc 2016;11:905-19. [CrossRef]
40. Morris GM, Huey R, Lindstrom W, Sanner MF, Belew RK, Goodsell DS, et al. AutoDock4 and AutoDockTools4:Automated docking with selective receptor flexibility. J Comput Chem 2009;30:2785-91. [CrossRef]
41. Collier TA, Piggot TJ, Allison JR. Molecular dynamics simulation of proteins. Methods Mol Biol 2020;2073:311-27. [CrossRef]
42. Tumskiy RS, Tumskaia AV. Multistep rational molecular design and combined docking for discovery of novel classes of inhibitors of SARS-CoV-2 main protease 3CLpro. Chem Phys Lett 2021;780:138894. [CrossRef]
43. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, et al. GROMACS:High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015;1-2:19-25. [CrossRef]
44. Berendsen HJ, van der Spoel D, van Drunen R. GROMACS:A message-passing parallel molecular dynamics implementation. Comput Phys Commun 1995;91:43-56. [CrossRef]
45. Huang W, Lin Z, van Gunsteren WF. Validation of the GROMOS 54A7 force field with respect to b-peptide folding. J Chem Theory Comput 2011;7:1237-43. [CrossRef]
46. Schüttelkopf AW, van Aalten DM. PRODRG:A tool for high-throughput crystallography of protein-ligand complexes. Acta Crystallogr D Biol Crystallogr 2004;60:1355-63. [CrossRef]
47. Lemkul JA. From proteins to perturbed hamiltonians:A suite of tutorials for the GROMACS-2018 molecular simulation package [article v1.0]. Living J Comput Mol Sci 2019;1:5068. [CrossRef]
48. Parrinello M, Rahman A. Polymorphic transitions in single crystals:A new molecular dynamics method. J Appl Phys 1981;52:7182-90. [CrossRef]
49. Oostenbrink C, Villa A, Mark AE, van Gunsteren WF. A biomolecular force field based on the free enthalpy of hydration and solvation:The GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem 2004;25:1656-76. [CrossRef]
50. Humphrey W, Dalke A, Schulten K. VMD:Visual molecular dynamics. J Mol Graph 1996;14:33-8, 27-8. [CrossRef]
51. Soman SS, Sivakumar KC, Sreekumar E. Molecular dynamics simulation studies and in vitro site directed mutagenesis of avian beta-defensin Apl_AvBD2. BMC Bioinformatics 2010;11 Suppl 1:S7. [CrossRef]
52. Kumari R, Kumar R, Lynn A, Open Source Drug Discovery Consortium. g_mmpbsa--a GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inf Model 2014;54:1951-62. [CrossRef]
53. Shaik NA, Hakeem KR, Banaganapalli B, Elango R. Essentials of Bioinformatics. Understanding Bioinformatics:Genes to Proteins. Vol. 1. Germany:Springer;2019. [CrossRef]
54. Daina A, Michielin O, Zoete V. SwissADME:A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci Rep 2017;7:42717. [CrossRef]
55. Lipinski CA, Lombardo F, Dominy BW, Feeney PJ. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliv Rev 2001;46:3-26. [CrossRef]
56. Veber DF, Johnson SR, Cheng HY, Smith BR, Ward KW, Kopple KD. Molecular properties that influence the oral bioavailability of drug candidates. J Med Chem 2002;45:2615-23. [CrossRef]
57. Perez A, Morrone JA, Simmerling C, Dill KA. Advances in free-energy-based simulations of protein folding and ligand binding. Curr Opin Struct Biol 2016;36:25-31. [CrossRef]
58. Raj S, Sasidharan S, Dubey VK, Saudagar P. Identification of lead molecules against potential drug target protein MAPK4 from L. donovani:An in-silico approach using docking, molecular dynamics and binding free energy calculation. PLoS One 2019;14:e0221331. [CrossRef]
59. Broni E, Kwofie SK, Asiedu SO, Miller WA 3rd, Wilson MD. A molecular modeling approach to identify potential antileishmanial compounds against the cell division cycle (cdc)-2-related kinase 12 (CRK12) receptor of Leishmania donovani. Biomolecules 2021;11:458. [CrossRef]
60. Gutiérrez-Rebolledo GA, Drier-Jonas S, Jiménez-Arellanes MA. Natural compounds and extracts from Mexican medicinal plants with anti-leishmaniasis activity:An update. Asian Pac J Trop Med 2017;10:1105-10. [CrossRef]
61. Passero LF, Brunelli ED, Sauini T, Pavani TF, Jesus JA, Rodrigues E. The potential of traditional knowledge to develop effective medicines for the treatment of leishmaniasis. Front Pharmacol 2021;12:690432. [CrossRef]
62. Azizi K, Shahidi-Hakak F, Asgari Q, Hatam GR, Fakoorziba MR, Miri R, et al. In vitro efficacy of ethanolic extract of Artemisia absinthium (Asteraceae) against Leishmania major L. using cell sensitivity and flow cytometry assays. J Parasit Dis 2016;40:735-40. [CrossRef]