Inhibiting Leishmania donovani Sterol Methyltransferase to Identify Lead Compounds Using Molecular Modelling

The recent outlook of leishmaniasis as a global public health concern coupled with the reportage of resistance and lack of efficacy of most antileishmanial drugs calls for a concerted effort to find new leads. The study combined In silico and in vitro approaches to identify novel potential synthetic small-molecule inhibitors targeting the Leishmania donovani sterol methyltransferase (LdSMT). The LdSMT enzyme in the ergosterol biosynthetic pathway is required for the parasite’s membrane fluidity, distribution of membrane proteins, and control of the cell cycle. The lack of LdSMT homologue in the human host and its conserved nature among all Leishmania parasites makes it a viable target for future antileishmanial drugs. Initially, six known inhibitors of LdSMT with IC50 < 10 μM were used to generate a pharmacophore model with a score of 0.9144 using LigandScout. The validated model was used to screen a synthetic library of 95,630 compounds obtained from InterBioScreen limited. Twenty compounds with pharmacophore fit scores above 50 were docked against the modelled three-dimensional structure of LdSMT using AutoDock Vina. Consequently, nine compounds with binding energies ranging from −7.5 to −8.7 kcal/mol were identified as potential hit molecules. Three compounds comprising STOCK6S-06707, STOCK6S-84928, and STOCK6S-65920 with respective binding energies of −8.7, −8.2, and −8.0 kcal/mol, lower than 22,26-azasterol (−7.6 kcal/mol), a known LdSMT inhibitor, were selected as plausible lead molecules. Molecular dynamics simulation studies and molecular mechanics Poisson–Boltzmann surface area calculations showed that the residues Asp25 and Trp208 were critical for ligand binding. The compounds were also predicted to have antileishmanial activity with reasonable pharmacological and toxicity profiles. When the antileishmanial activity of the three hits was evaluated in vitro against the promastigotes of L. donovani, mean half-maximal inhibitory concentrations (IC50) of 21.9 ± 1.5 μM (STOCK6S-06707), 23.5 ± 1.1 μM (STOCK6S-84928), and 118.3 ± 5.8 μM (STOCK6S-65920) were obtained. Furthermore, STOCK6S-84928 and STOCK6S-65920 inhibited the growth of Trypanosoma brucei, with IC50 of 14.3 ± 2.0 μM and 18.1 ± 1.4 μM, respectively. The identified compounds could be optimised to develop potent antileishmanial therapeutic agents.


Introduction
The challenges associated with control coupled with the frequent epidemic outbreak have rendered leishmaniasis an important public health problem [1][2][3]. Consequently, the presence of the ugly lesions after treatment and the severity of those caused by Leishmania donovani and Leishmania infantum in humans has rendered visceral leishmaniasis a global public health concern that needs urgent attention [1,3]. Currently, there are no vaccines and the number of chemotherapeutic options is limited [4,5]. Drugs including pentamidine, amphotericin B, paromomycin, miltefosine, and stibogluconate are ineffective and toxic, and the parasite has developed resistance, requiring renewed and concerted efforts to identify and develop new antileishmanial chemotypes [6,7]. Some efforts have been made towards finding new treatment options from natural products [8,9]; however, their structures and chemical synthesis are complex [10,11].
Synthetic small molecules with drug-like properties have been explored for the treatment of various ailments including leishmaniasis [12,13]. Recent studies show that about 60% of the drugs used for treating most diseases originated from synthetic small molecules [14,15]. The inherent properties of synthetic small molecules such as the ability to cross biological barriers and to modulate diverse biological targets render them viable drug candidates [16].
Target identification and validation are pivotal in drug design and the success of drug development depends primarily on choosing the right target. Studies in the design of antileishmanial compounds have identified some novel pathways and protein targets necessary for the parasite's survival [17,18]. Some of these targets have been validated [19,20], while investigation on others is still ongoing. Ergosterol biosynthesis is involved in various biological functions including plasma membrane formation, membrane fluidity, distribution of membrane proteins, and control of the cell cycle [21]. Ergosterol is essential for optimal mitochondrion function in Leishmania parasites [22]. This pathway is catalysed by several enzymes with homologues in the human host and hence has the potential to cause off-target effects, mainly as a result of poor selectivity [23]. Sterol methyltransferase, which catalyses the transfer of the methyl group from S-adenosine-methionine to the C24 position of the sterol side chain during ergosterol biosynthesis, lacks a homologue in humans [22,24], and is thus considered an essential drug target.
Reports have shown that 22,26-azasterol suppresses the growth of L. donovani with a half-maximal inhibitory concentration (IC 50 ) of 8.9 µM [25]. To improve the inhibitory effect of 22,26-azsterol, some analogues ( Figure 1) were generated, and their activities were found to be dependent on the presence of ammonium or sulfonium functionality in the side chain [25]. Similarly, ezetimibe, imipramine, and simeprevir ( Figure 1) have shown promising inhibitory activity against L. donovani by making morphological changes to the plasma and mitochondrion membranes with IC 50 of 30, 28.6, and 51.49 µM, respectively [26][27][28]. A fragment-based de novo design was used to predict potential inhibitors against L. donovani sterol methyltransferase (LdSMT) [29]. Though several hit compounds have been identified targeting this enzyme, only a few progressed to the clinical evaluation stage [30][31][32][33][34].
In silico techniques offer alternative platforms for the identification of novel hits in a timely and cost-effective manner [35,36]. These approaches have improved rational drug design and augmented the identification of potentially new drug candidates. This study, therefore, employs ligand-based pharmacophore virtual screening, molecular docking, molecular dynamics (MD) simulations, molecular mechanics Poisson-Boltzmann surface area (MM/PBSA), biological activity predictions, and in vitro studies to support the identification of potential inhibitors. The synthetic compound library was screened against LdSMT followed by the elucidation of the mechanisms of binding of the proposed molecules.

Materials and Methods
The schematic workflow employed in the project is shown in Figure 2. First, six known inhibitors of LdSMT with IC 50 less than 10 µM were used to generate a pharmacophore model employing the merged feature embedded in LigandScout version 4.3 [37]. The pharmacophore model was then used to screen a library of 95,630 synthetic compounds from InterBioScreen limited. Ligands with pharmacophore fit scores above 50 were docked against a previously modelled structure of LdSMT [29,38]. The selected hit compounds were evaluated via physicochemical, pharmacological, and toxicity profiles, as well as biological activity predictions, MD simulations, and MM/PBSA computations. Three potential lead compounds were then screened in vitro to validate the In silico studies.

Target Preparation
The three-dimensional (3D) structure of LdSMT previously elucidated using Modeller version 10.2 [29,38] was used. Energy minimisation of the target protein was carried out by employing the Optimized Potential for Liquid Simulations (OPLS)/All Atom force field [39] via the Groningen Machine for Chemical Simulations GROMACS version 2018 (GROMACS 2018) [40,41]. Biovia Discovery Studio Visualizer v.19.1.0.18287 [42] was used to visualise the energy-minimised structure, remove water molecules, and solvate the protein, and then the result was saved in Protein Data Bank (pdb) format, which was later converted to AutoDock Vina's [43] compatible pdbqt format.

Binding Site Prediction
The amino acid residues and the probable volume and area of the binding site of the modelled protein were determined using the Computed Atlas of Surface Topography of proteins (CASTp) [44] and Biovia Discovery Studio Visualizer v.19.1.0.18287 [42], as indicated in a previous study [29].

Ligand-Based Pharmacophore Virtual Screening
Six of the known inhibitors ( Figure 1) of LdSMT with IC 50 less than 10 µM were employed in generating the pharmacophore model using LigandScout version 4.3 [37]. The structure data file (sdf) formats of the ligands were used via LigandScout's Ligand-Based Modeling Perspective v.4.3 [37]. The default settings of OMEGA best were used to generate a maximum of 200 conformations per ligand.

Retrieval and Preparation of Chemical Library for Pharmacophore-Based Screening
A chemical library of 95,630 synthetic compounds was retrieved from InterBioScreen limited [45] and used for the pharmacophore-based virtual screening.

Pharmacophore-Based Virtual Screening of the Libraries
The 95,630 chemical entities retrieved from InterBioScreen limited were utilised for the pharmacophore screening using LigandScout v.4.3 [46], first by converting from "sdf" to "lbd" and then by screening against a validated pharmacophore model.

Validation of Pharmacophore Model and AutoDock Vina
Before performing the virtual screening, both the pharmacophore model and AutoDock Vina v.1.2.0 [43] were validated to assess their potential for distinguishing between active compounds and decoys.

Pharmacophore Model Validation
The receiver operating characteristic (ROC) curve and enrichment factors (EFs) were used in validating the pharmacophore model in LigandScout v.4.3 [37]. The six inhibitors of SMT and their decoys were used to generate the ROC curve from where the area under the curve (AUC) and EFs were calculated.

Validation of AutoDock Vina
The simplified molecular input line entry system (SMILES) of six of the known inhibitors with IC 50 less than 10 µM served as active compounds for the generation of 50 decoys each using the Directory of Useful Decoys (DUD-E) [47]. The ROC curves generated via EasyROC v.1.3.3 [48] were used to validate the docking protocol of AutoDock Vina v.1.2.0 [43].

Molecular Docking Studies of Chemical Entities with Good Pharmacophore Fit Scores
AutoDock Vina [43] interfaced with PyRx v.0.8 [49], which employs empirical scoring functions, was used for the virtual screening. Ligands with good pharmacophore fit scores were obtained after screening the pharmacophore model against the library combined with 22,26-azasterol and the three known drugs (amphotericin B, miltefosine, and paromomycin). The energy-minimised protein (LdSMT) was used for molecular docking via AutoDock Vina v.1.2.0 [43]. The charge, hydrogen bond network, and histidine proto-nation state of the protein were assigned after pdbqt conversion [32]. The grid box size (91.445 × 73.502 × 78.352) Å 3 with the centre (72.200, 58.009, 13.302) Å was used with an exhaustiveness of 8 [29].

Characterisation of Mechanism of Binding
Biovia Discovery Studio Visualizer v.19.1.0.18287 [42] was employed in the elucidation of the protein-ligand interactions.

ADMET Properties and Drug-Likeness Assessment
The absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties were determined using SwissADME [50] and the OSIRIS Property Explorer in DataWarrior [51]. This was carried out to evaluate the drug-likeness and the pharmacologic properties of the selected compounds.

Prediction of Biological Activity of Selected Compounds
The biological activity of the compounds was predicted using Prediction of Activity Spectra for Substance (PASS) [52]. The SMILES files of the molecules were used as inputs.

Molecular Dynamics Simulation Analyses of Protein and Protein-Ligand Complexes
A 100 ns MD simulation was adopted [29] for LdSMT and LdSMT-ligand complexes via a sample interval for configuration of 100 ps using GROMACS 2018 [40,41]. Initially, the topology and position restrain files for the protein were generated using OPLS/All Atom force field [39]. LigParGen [53] was employed to generate the topology and parameter files of the ligands. The unbound protein and protein-ligand complexes were centred in a cubic simulation box with a 1 nm distance from the edges to restrain particle movements that have the potential to cause effects on the surface atoms [54]. System solvation was carried out using a single-point charge and then neutralised to achieve a molarity of 0.15 M using NaCl. Each system was energy-minimised for 50,000 steps using the steepest descent method. A 100 ps equilibration was carried out using several substances, volume and temperature, and several substances, pressure, and temperature ensembles [54] to ensure a well-equilibrated system was achieved at 300 K and a pressure of 1 bar. The root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), the radius of gyration (Rg), and surface accessible surface area (SASA) trajectories were plotted using Qtgrace [55]. MM/PBSA computations of the complexes were carried out using g_MM/PBSA, which calculates the free binding energy and the individual energy contributions of the residues [56]. were used for the cell viability assay. The promastigotes were grown at 25 • C in M199 medium supplemented with 10% heat-inactivated foetal bovine serum, penicillin G sodium (100 g/mL), and streptomycin sulphate (100 g/mL), and subcultured every 72 h in the same medium at a density of 2 × 10 5 cells/mL.
To determine the concentration of a compound that inhibits the growth of 50% (IC 50 ) of the Leishmania parasite population, promastigotes at a density of 2 × 10 5 cells were incubated in triplicate without or with the compounds at varying concentrations (100-0.781) µg/mL and kept at 25 • C for 72 h. The Alamar blue test was used to determine the vitality of the parasites as previously described [57]. At an excitation wavelength of 530 nm and an emission wavelength of 590 nm, fluorescence was measured using a Varioskan Lux Elisa microplate reader (Thermo Fisher Scientific, Waltham, MA, USA). The fluorescence counts were plotted against drug concentrations, and the IC 50 was calculated using a dose-response curve via GraphPad Prism version 9 (GraphPad Software Inc., San Diego, CA, USA). The reference drug used in the assay is amphotericin B, a known drug for treating leishmaniasis.

In Vitro Culture and Cell Viability Assay for Trypanosoma brucei brucei GUTat 3.1
Similar to the approach described in Section 2.12, to determine the IC 50 of the compounds against Trypanosoma brucei brucei GUTat 3.1 parasite population, the modifications were trypanosomes at a density of 4×10 3 cells, which were kept at 37 • C for 72 h [58]. The control drug used for this assay was diminazene aceturate (Sigma-Aldrich, Kent, UK), a veterinary drug used to combat infections of trypanosomes in cattle.

Results and Discussion
The results from the study are presented in the following sections, including the pharmacophore-based design, molecular docking, ADMET profiling, prediction of biological activities of selected molecules, MD simulations, MM/PBSA computations, and biological activity evaluations.

Predictions of LdSMT Binding Site Residues
A binding site is a pocket or groove on a protein that serves as a site for a ligand or biological macromolecule to bind with specificity. Ligands bind at active sites to induce biochemical reactions including methylation, demethylation, and phosphorylation. CASTp [44], which employs Delaunay triangulation, alpha shape, and discrete flow methods to identify topographic features and measure areas and volumes, was used to predict the binding site of LdSMT in previous work [29]. The predicted binding site volume was 446.632 Å 3 and the area was 905.262 Å 2 , and the critical amino acid residues elucidated included Ala30, Asp108, Val109, Tyr206, Gly207, Trp208, Met210, Tyr220, Tyr275, Leu278, Glu324, Ser328, Leu329, Val330, and Val331 [29].

Pharmacophore Generation
For pharmacophore generation and the prediction of features, six LdSMT inhibitors (22,26-azasterol and X1-X5) ( Figure 1) with IC 50 < 10 µM were used. In addition, all six ligands and their decoys were employed to validate the pharmacophore model using the merged feature of LigandScout version 4.3 [37], similar to previous studies [59]. A predicted score of 0.9144 was generated with features consisting of three hydrophobic interactions (H), one positive ionisable (PI), one hydrogen bond acceptor (HBA), and one hydrogen bond donor (HBD) (Figure 3).

Validation of Generated Pharmacophore Model
The ROC curve evaluates the ability of the model to efficiently distinguish between a compendium of active and inactive compounds [59]. The ROC curve was obtained by plotting a graph of 1-specificity (Equation (1)) on the x-axis against sensitivity (obtained using Equation (2)) on the y-axis [60]. Notably, 1-specificity is called false positive fraction (FPF) and the sensitivity is also known as the true positive fraction (TPF).
where FP denotes false positive, TN as true negative, TP as true positive, and FN as false negative. Consequently, while a poor test has a ROC curve falling on the diagonal line, a perfect test has the ROC curve passing through the upper left corner with approximately 100% sensitivity and 100% specificity [60]. The ROC curve (Figure 4a) falls to the upper left corner, implying a reasonably valid pharmacophore model. The AUC ranges from 0 to 1, with values close to 1 signifying a perfect model with the ability to distinguish active from inactive compounds [59]. The six known inhibitors serving as active compounds were used to generate 300 decoys made up of 50 decoys from each. A total of 306 compounds comprising six active molecules and 300 decoys were screened via the 3D pharmacophore model and the AUC was computed as 1 ( Figure 4a). The pharmacophore-guided approach was used to identify inhibitors of Mycobacterium ulcerans Cystathione γ-synthase MetB and an AUC of 0.7 was obtained after validation [59]. In addition, EF, which is the ratio of the number of identified active compounds within a pool of highly ranked hits to those randomly selected from the original dataset, was computed to evaluate the effectiveness of the pharmacophore model [61,62]. The EF was computed using Equation (3) [61]. EF = (N hits sampled/ N sampled )/(N total hits/ N total ) (3) where N hits sampled is the number of identified true hits in the hit list, N sampled is the number of compounds in the hit list, N total hits is the number of hits in the dataset, and N total is the total number of compounds in the dataset. A pharmacophore model with EF less than 1 suggests the classification ability is low. Similarly, EF equal to one shows that the model has equal chances of classifying active and inactive compounds. In addition, an EF greater than one suggests the model has a greater potential of classifying active compounds. A pharmacophore model generated for the design of sigma-1 ligands was classified as poor due to the computed EF being close to 1 at 10% [63]. Similarly, a model with EF score above 1 classified the inhibitors significantly better [61]. Consequently, the EFs were obtained as 52.0, 26.0, 26.0, and 26.0 for 1%, 5%, 10%, and 100%, respectively, implying the effectiveness of the pharmacophore model in classifying hits amongst the dataset.

Validation of Docking Protocol
Most docking algorithms, including the Lamarckian genetic algorithm and the empirical scoring for free energy binding employed in AutoDock Vina, sometimes fail to accurately predict the pose and scoring functions, warranting validation [64]. The com-puted AUC under the ROC was used to evaluate the performance of AutoDock Vina v.1.2.0 [43] in distinguishing the active from inactive compounds. To validate the docking protocol, the six active compounds and 300 decoys were virtually screened against the LdSMT to generate the ROC curve and the AUC was computed as 0.9997 (Figure 4b) [65,66]. The ROC curve (Figure 4b) aligns with the upper left corner of the graph, which signifies a reasonably excellent docking protocol. The AUC of 0.9997 suggests that AutoDock Vina v.1.2.0 [43] has a plausible ability to distinguish active compounds from decoys.

Pharmacophore-Based Virtual Screening of Library
Pharmacophore-based virtual screening is the generation of a 3D query by exploring the chemical features responsible for the biological activity from a reasonable number of structures for the identification of new chemotypes [67,68]. The surge in interest is due to the high enrichment of novel actives generated from a library of chemical entities [67]. The validated pharmacophore model was used as a 3D query to screen a chemical library of 95,630 synthetic compounds from InterBioScreen limited with a pharmacophore fit score of 50 as the threshold. STOCK6S-64914 and STOCK6S-47366 had the highest and least pharmacophore fit scores of 59.39 and 57.15, respectively, among the 20 shortlisted compounds (Table 1).

Molecular Docking Analysis
Molecular docking predicts the orientation and conformation of the ligand within the binding site of the target protein. Strategically, molecular docking is supposed to mimic a biological system, but being a computational-based approach [69], it often suffers certain drawbacks. In the current study, the effects and implications of solvent present in the protein that otherwise may affect the pose of the ligand [70] were not considered, as water molecules were deleted during protein preparation. MD simulation studies were carried out to ascertain the contribution of the solvent in the free binding energy of the protein-ligand complexes similar to previous studies [71,72]. In addition, the preparations of both the protein and the ligand before docking have the potential to introduce errors, leading to missing bonds and abnormal geometries. This occurs when converting from pdb or sdf to the AutoDock Vina format pdbqt. Visualisation analysis of the protein and ligands after preparation and file conversion before docking ensures that this does not affect the pose of the ligands. The exhaustiveness in the current study was set to default 8, which statistically may not represent the actual number of possible conformations of ligands in the binding pocket of the biological target. This is normally occasioned to save computational time and energy. Furthermore, the top-ranked pose was chosen as the best pose with the least binding energy. Interestingly, this is not always the case, as the top-ranked score may not necessarily represent the best pose. Despite these drawbacks, molecular docking has provided an ideal platform for studying molecular inclusions at the atomic level [73]. Molecular docking has been explored in the identification of lead compounds including gentisic acid 5-O glucoside [59], simeprevir [28], and bisindolylmaleimide IX [74] by targeting Buruli ulcer [59], visceral leishmaniasis [28], and SARS-CoV-2 [74], respectively. In all these instances, the binding energies associated with the protein-ligand complexes were used as part of the criteria to shortlist the hits. The 20 shortlisted compounds from the pharmacophore screening were docked into the catalytic domain [29] of the LdSMT structure with binding energies comparable to or lower than that of 22,26-azasterol (−7.6 kcal/mol). The binding energies of nine compounds (Table 2) were also lower than those of the antileishmanial drugs amphotericin B (−5.3 kcal/mol), miltefosine (−5.0 kcal/mol), and paromomycin (−4 kcal/mol) ( Table 2). The binding energies ranged from −7.5 to −8.7 kcal/mol with STOCK6S-07353, STOCK6S-16994, and STOCK6S-06707 having −7.5, −7.5, and −8.7 kcal/mol, respectively. The results showed that the nine selected compounds possess better binding energies and can be explored as further potential inhibitors.

Protein-Ligand Interaction
Molecular recognition due to protein-ligand interactions is the basis of most cellular mechanisms. The examination of this at the atomistic level provides insights into the design of small-molecule drugs. Molecular docking was undertaken to elucidate the mechanism of binding of the ligands as well as to identify small molecules with high cooperativity [75] and affinity to the LdSMT. An evaluation of the 2D interactions [42] showed that the compounds including 22,26-azasterol and the three drugs interacted with amino acid residues Asp25, Ala28, Phe93, Phe100, Trp208, Asp218, Arg227, Ile228, Ala311, Glu320, and Leu322, consistent with previous studies [29,76]. Specific interactions in all of the complexes included pi-anion, pi-pi stacking, pi-alkyl, pi-sigma, carbon-hydrogen, and hydrogen bonds.

ADMET Profiling
The drug-likeness of the nine selected compounds was evaluated using Lipinski's rule of five (Ro5) [78,79]. According to the Ro5, a compound is classified as drug-like and orally active if it conforms to the following criteria: a molecular weight ≤500 Da, an octanol-water partition coefficient log p ≤ 5, HBD ≤ 5, and HBA ≤ 10 [79]. The Ro5 has been applied in many studies, especially during the shortlisting of hits [80][81][82]. The physicochemical properties and drug-likeness computed via SwissADME [50] showed that none of the compounds violated any of the Ro5 (Table 3). Orally active drugs, according to Veber's rule, must possess a maximum of 10 rotatable bonds as well as a polar surface area of less than 140 Å 2 [83]. All nine hits were predicted to have less than 10 and 140 Å 2 for the number of rotatable bonds and polar surface area, respectively. STOCK6S-64941 had the highest number of rotatable bonds (6) and a polar surface area of 80.57 Å 2 . The results (Table 3) suggest that the selected compounds possess the potential to be orally active [83].
Molar refractivity (MR) measurement was conducted to assess the movement of the compounds through the body and their ability to reach their targets of inhibition at optimum concentrations [84]. MR for druggable candidates is recommended to be between 40 and 130 [29]. All of the selected compounds were found to possess MR within this range except STOCK6S-64941 with an MR of 130.16 (Table 3).
For drugs to be transported by the circulatory system to their site of activity, they are required to be soluble [85]. The computed solubility profile of the compounds showed them as soluble and hence possess the potential to reach their target sites at maximum concentrations (Table 3). Furthermore, the bioavailability score (BS), which measures the ability of drugs to reach systemic circulation, was also estimated [86]. The computed BS for the compounds was 0.55, similar to the reported studies (Table 3) [29,86]. Table 3. Predicted physicochemical and pharmacological properties of the selected hit compounds. Molecular weight (MW), number of rotatable bonds (NRB), molar refractivity (MR), topological polar surface area (TPSA), gastrointestinal absorption (GI), blood-brain barrier (BBB), number of Ro5 violations (vRoF), bioavailability score (BS), and solubility score (SC) are presented. The ability of a drug to cross the blood-brain barrier (BBB) to the brain is called BBB permeation [87]. Drugs with the ability to cross the BBB could bind to certain receptors to elicit the required pharmacological activities within the brain parenchyma. STOCK6S-55084, STOCK6S-64941, and STOCK6S-19430 (Table 3) were predicted not to be BBB permeants. Previous reports showed that Leishmania parasites infect and inflame the central nervous system (CNS) [88][89][90] and drugs with the potential of crossing the BBB could mediate in the elimination of the parasites in the CNS.

Compound
Next, we investigated the gastrointestinal absorption of the compounds. Together with 22,26-azatserol, they were predicted to possess high gastrointestinal absorption, suggesting a high probability of absorption into the bloodstream. On the contrary, the three drugs, amphotericin B, miltefosine, and paromomycin, were predicted to show low gastrointestinal absorption (Table 3). A previous study to identify potential inhibitors of L. donovani cell division cycle-2-related kinase 12 (CRK12) predicted compounds to possess high gastrointestinal absorption [32].
Furthermore, the toxicity profile (Table 4) of the compounds was evaluated using OSIRIS DataWarrior 5.0.0 [51]. The results showed that all of them were predicted not to be mutagenic or tumorigenic. On the other hand, STOCK6S-0670 was predicted as a high irritant, whilst STOCK6S-65920 and STOCK6S-55084 were found to possess low reproductive effects. Compounds STOCK6S-19430, STOCK6S-14893, and STOCK6S-16994 were predicted to possess high irritant and reproductive effects (Table 4).

Biological Activity Predictions of Selected Hit Compounds
The biological activities of the selected compounds were predicted [52] based on the structure-activity relationship. Compounds predicted with a probability of activity (Pa) greater than the probability of inactivity (Pi) have high chances of exhibiting the expected biological activities, hence could be considered for experimental validation [52,81].
STOCK6S-64941 and STOCK6S-07353 were predicted as potential aspulvinone dimethylallyl transferase and indole pyruvate C-methyltransferase inhibitors (Supplementary Table S1), respectively. The LdSMT belongs to the family of methyltransferases to which aspulvinone dimethylallyl transferase and indole pyruvate C-methyltransferase also belong [95][96][97]. The predicted inhibition of these two targets suggests the potential of the selected hit compounds to suppress L. donovani sterol methyltransferase.
Furthermore, five of the selected hits comprising STOCK6S-06707, STOCK6S-55084, STOCK6S-19430, STOCK6S-14893, and STOCK6S-16994 were predicted for Alzheimer's treatment (Supplementary Table S1) with Pa greater than Pi. A recent in vitro study identified perphenazine, a known drug for Alzheimer's treatment, to possess antileishmanial potential with EC 50 of 1.2 µg/mL [98]. Similarly, clomipramine, a drug for Alzheimer's treatments, induced programmed cell death in L. amazonensis via a mitochondrial pathway disruption with an IC 50 of 8.31 µM [99,100]. The predicted biological activities and mechanisms suggest that the selected compounds have potential antileishmanial scaffolds.

MD Simulations Analyses
Three of the hit compounds, STOCK6S-06707, STOCK6S-84928, and STOCK6S-65920, were selected for downstream MD simulation analysis. They possessed low binding energies, reasonable pharmacological profiles, and desirable biological activities. The control used for the MD was 22,26-azasterol. A 100 ns MD simulation was performed to assess the structural stability and the conformational dynamics of the LdSMT-analogue complexes. The RMSDs were computed to ascertain the stability of the protein-ligand complexes (Supplementary Figure S2). The complexes showed comparable stability to the unbound protein and protein-22,26-azasterol complex. The unbound LdSMT and the protein-ligand complexes experienced a high volatility within the first 10 ns of the MD simulation with an RMSD fluctuation of 0.15 nm (Supplementary Figure S2a), which is classified as acceptable, since the drifting falls within the range 0.1-0.3 ns [101]. Equilibration was therefore reached afterward with a low RMSD fluctuation of ≤0.05nm. Remarkably, the RMSD results produced were corroborated by the RMSF values (Supplementary Figure S3a).
The rigidness of the unbound protein and the complexes for the entire simulation period is determined by Rg [32]. A folded protein normally maintains a relatively steady Rg over a given simulation period when compared to the unfolded period. The computed Rg plot showed that the complexes maintained steady stability over the entire period with the LdSMT-STOCK65-06707 and LdSMT-STOCK6S-84928 showing the lowest average Rg of 1.975 ± 0.02 and 1.975 ± 0.01 nm, respectively (Supplementary Figure S4a). Similarly, 22,26-azasterol, STOCK6S-65920 complexes, and the unbound LdSMT recorded respective averages of Rg of about 2.0 ± 0.01, 1.975 ± 0.03, and 2.0 ± 0.01 nm, respectively. As per the PDF analysis, all of the complexes and the native protein had only one peak. Both the unbound LdSMT and LdSMT-ligand complexes had an average Rg of 2 nm (Supplementary Figure S4b). Furthermore, the Rg graph ( Figure 6) indicates that the compactness of the LdSMT-STOCK6S-06707, LdSMT-STOCK6S-84928, and LdSMT-STOCK6S-65920 was maintained and kept steady after complex formation.
To determine the effect of loose packing on solvent accessibility, the SASA of the entire protein was computed to provide information on the protein-solvent interactions [103,105]. Generally, high SASA values imply an increase in structural enlargement for the proteinligand complexes under the influence of solvent surface charges, resulting in a more flexible and unstable conformation [105]. A complex of L. donovani 3-mercapto pyruvate sulfurtransferase and rutin complex had an average SASA value of 44 nm 2 , which was more compact than the unbound protein with an average SASA of 48 nm 2 [103]. In the current study, no differences were observed (Supplementary Figure S5a  To better understand the binding interactions that influenced the binding energy of the various protein-ligand complexes, the MD simulation was carried out to calculate the number of HB formed for all the LdSMT-analogue complexes. The number of HB formed were (6 ± 1), (5 ± 1), (3 ± 1), and (2 ± 1) ( Figure 6) for the LdSMT-STOCK6S-06707, LdSMT-22,26azasterol, LdSMT-STOCK6S-84928, and LdSMT-STOCK6S-65920 complexes, respectively. Altogether, the hydrogen and hydrophobic interactions showed strong binding of the hit compounds in the active site of the LdSMT.

MM/PBSA Computation of Free Binding Energies
Comparatively, binding free energy calculations from simulation studies have been proven to be more accurate than their docking counterparts and can prioritise compounds more efficiently for experimental evaluation [106]. Therefore, this work further employed MM/PBSA in computing the binding free energies of the protein-ligand complexes. Notably, a high negative binding free energy of a complex signifies a strong affinity of the ligand to the target protein. Among the selected hit compounds, STOCK6S-06707 with the lowest binding energy of −8.7 kcal/mol from the docking studies also recorded the lowest binding free energy of −371.146 ± 2.105 kJ/mol (Table 5). While STOCK6S-84928 and STOCK6S-65920 showed binding free energies of −129.725 ± 4.799 and −149.899 ± 3.6 kJ/mol, respectively, their binding free energies were higher than 22,26azasterol (−221.527 ± 3.716 kJ/mol). The lead compounds showed promising binding free energies comparable to 22,26-azasterol and therefore warrant experimental validation. The free binding energy was contributed by van der Waal (∆G vdW ), electrostatic (∆G ele ), polar solvation (∆G pol, sol ), and SASA (∆G SASA ). Among these energies, van der Waal, responsible for the embedded hydrophobic interactions between the residues in the binding pocket and lead compounds, was the major contributor to the binding free energy (Table 5), followed by electrostatic. The results presented in this work corroborate previous studies that have shown that both van der Waal and electrostatic energies were responsible for the stability of protein-ligand complexes [29,32]. Despite the role of entropy in computing the free energy of binding to protein-ligand complexes [107], the current study did not take into account this contribution, because entropy tends to produce negligible effects when comparing the binding strength of ligands inhibiting the same protein target [108] and could be computationally expensive. Favourable interactions between the binding site residues of the target protein and ligand result in the stability of the complex, and this was assessed by the per-residue decomposition, which computes the energy contribution of each of the amino acids. Generally, amino acids that contribute energies greater than +5 kJ/mol or lower than −5 kJ/mol are regarded as critical for ligand binding [59]. From the docking studies, residues Asp25, Ala28, Phe93, Phe100, Trp208, Asp218, Arg227, Ile228, Ala311, Glu320, and Leu322 were found to be key for ligand binding. The MM/PBSA per-residue decomposition computations revealed that Asp25 and Trp208 (Figure 7) contributed energies of −20 and −48 kJ/mol, respectively, for the stability of the LdSMT-STOCK6S-06707 complex. Furthermore, the amino acid Trp208 contributed an energy of −10 kJ/mol () for the stability of STOCK6S-65920. The 22,26-azasterol complex generated per-residue decomposition trajectory with Glu102, Phe194, Phe198, and Gly200 (Supplementary Figure S6b,c), contributing energies above and below the threshold critical for ligand binding. Similarly, the per-residue decomposition trajectory of the STOCK6S-84928 complex shows Trp208 to produce energy of −7 kJ/mol for ligand binding. An earlier study to identify the potential inhibitors of LdSMT showed similar residues as being critical for ligand binding [29]. The amino acid residues within the binding pocket, including Asp25 and Trp208 of LdSMT, contributed meaningful energies to the stability of the complexes.

In Vitro Evaluation of Identified Hits against Leishmania donovani
The antileishmanial properties of three purchasable hit compounds, STOCK6S-06707, STOCK6S-84928, and STOCK6S-65920, were tested in vitro. The effects of different concentrations (Supplementary Figure S7) on the survival of L. donovani promastigotes were studied with amphotericin B treatment as a positive control (Table 6). After 72 h of treatments of L. donovani promastigote cultures with varying concentrations of the three compounds, effect-based dose finding analysis revealed that concentrations of 21.9 ± 1.5 µM, 23.5 ± 1.1 µM, and 118.3 ± 5.8 µM (Table 6) for STOCK6S-06707, STOCK6S-65920, and STOCK6S-84928, respectively, eliminated 50 % of L. donovani promastigotes. Amphotericin B exhibited IC 50 of 6.56 ± 0.06 µM. Despite the significant inhibitory effects of the hit compounds, their potencies were lower compared to amphotericin B. This notwithstanding, the compounds can be optimised to improve potency.
Both Leishmania and Trypanosoma spp. belong to the family of kinetoplastids [109] and some compounds have been shown to inhibit the survival of both Leishmania parasites and Trypanosoma spp. [110][111][112]. Compound SQ10, which suppressed the growth of L. donovani, was also reported to possess activity against T. cruzi [112]. Moreover, 22,26azasterol, which inhibits sterol methyltransferase, also inhibits T. cruzi and T. brucei [25]. The anti-trypanosomal effects of the three compounds against T. brucei with diminazene as a positive control (Table 6)

Contribution to the Field
In addition to the limited number, the available chemotherapeutic agents for the treatment of visceral leishmaniasis are plagued with inefficiencies and chemoresistance [113]. The integrated In silico and in vitro study, therefore, augments the current efforts in discovering antileishmanial agents by identifying inhibitors with the potential to circumvent the activities of Leishmania donovani. The elucidation of the mechanisms of binding of the potential lead compounds provides the platform for the identification of other potent inhibitors for leishmaniasis treatment targeting sterol methyltransferase.

Conclusions
Sterol methyltransferase is expressed in all Leishmania parasites, but does not have human homologues, making it a potential target for drug design against leishmaniasis. This study used In silico techniques including pharmacophore-based virtual screening, molecular docking, and MD simulations, as well as in vitro approaches to identify potent inhibitors against LdSMT. Three compounds, STOCK6S-06707, STOCK6S-84928, and STOCK6S-65920, had a high binding affinity with binding energies lower than those of known inhibitors including 22,26-azasterol (−7.6 kcal/mol) and the antileishmanial drugs amphotericin B (−5.3 kcal/mol), miltefosine (−5.0 kcal/mol), and paromomycin (−4.0 kcal/mol). Molecular dynamics simulations and mechanistic studies revealed the leads to interact with critical amino acids in the binding site in the attenuation of LdSMT. The compounds were predicted as antileishmanial with acceptable pharmacological and toxicity profiles. STOCK6S-06707, STOCK6S-84928, and STOCK6S-65920 were shown to suppress promastigotes of L. donovani cultures with IC 50 values of 21.9, 23.5, and 118.3 µM, respectively. When the hits were tested for their anti-trypanosomal activity against T. brucei, STOCK6S-84928 and STOCK6S-65920 were active with IC 50 of 14.3 and 18.1 µM, respectively. The promising antileishmanial compounds identified could serve as a basis for the design of potent biotherapeutic moieties.