The protein sequence of LdSMT (SCMT1/GenBank ID: AAR92099.1) was retrieved from the National Center for Biotechnology Information (NCBI) database (Boratyn et al., 2013). To identify similar proteins as templates for building the 3D structure of LdSMT, the fast-all (FASTA) format of the sequence was aligned with the homologous sequences of crystal protein structures in the Protein Data Bank (PDB) (Berman et al., 2000) using the Basic Local Alignment and Sequencing Tool (BLAST) (Burley et al., 2017).

Due to the unavailability of an experimentally elucidated 3D structure for LdSMT, molecular modeling techniques were used to predict a reasonably accurate protein structure (Crentsil et al., 2020). Modeller version 10.2 (Eswar et al., 2008) was employed in building the 3D structure using three different templates: i) the 3BUS template was used to generate five models of the LdSMT (Eswar et al., 2008), ii) the protein structure with PDB ID 4PNE was also used to model five different structures of the LdSMT (Eswar et al., 2008), and iii) multitemplate homology modeling was employed to generate five structures of the LdSMT (Eswar et al., 2008). A total of three templates comprising 4PNE, 3BUS, and 6UAK were used in the multitemplate homology modeling approach. Modeller 10.2 was employed in generating all the potential structures of the LdSMT (Fiser and Šali, 2003; Eswar et al., 2008). The reasonably best model in each approach was selected based on the discrete optimized protein energy (DOPE) scores (Eswar et al., 2008).

The structural quality and accuracy of the best models from each approach were assessed using PROCHECK (Laskowski et al., 2012) with results reported as Ramachandran plots (Anderson et al., 2005). Further validation with VERIFY 3D (Sasin and Bujnicki, 2004), ERRAT (Colovos and Yeates, 1993; Messaoudi et al., 2013), and PROVE (Sasin and Bujnicki, 2004) was also performed. The reasonably best model was selected based on all the quality assessments for the molecular docking studies. Protein Structure Analysis (ProSA) (Wiederstein and Sippl, 2007) was then used to investigate the problematic regions of the selected model.

The plausible binding sites of the selected protein were determined using the Computed Atlas of Surface Topology of proteins (CASTp) (Dundas et al., 2006; Tian et al., 2018). The predicted sites were visualized using PyMOL (PyMOL Molecular Graphics System, Version 1.5.0.4, Schrödinger, LLC, New York, USA) (Lighthall et al., 2010) and Chimera 1.16 (Pettersen et al., 2004). The predicted binding sites with relatively small areas and volumes, where no ligand could fit, were ignored (Agyapong et al., 2021).

Shape similarity searching and pharmacophore screening were undertaken using the spatial data file (sdf) format of 22,26-azasterol via BRUSELAS (Banegas-Luna et al., 2019). A total of 100 ligands were generated with varying degrees of similarities to the input ligand, 22,26-azasterol. All ligands devoid of the steroidal core with varying degrees of similarities were selected for binding affinity prediction using AutoDock Vina v1.2.0 (Trott and Olson, 2010).

AutoDock Vina (Trott and Olson, 2010) was used for molecular docking. The molecular docking process was performed in two different stages. The first docking stage involved ligands obtained from the scaffold hopping process (Pathania et al., 2021), while the second involved compounds obtained from the de-novo design studies (Agyapong et al., 2021). Altogether, 1,448 ligands were used for the docking studies.

For the first stage, the ligands were obtained from the scaffold hopping and their structural derivatives fetched from the DrugBank (Wishart et al., 2018), PubChem (Kim et al., 2021), ZINC15 (Sterling and Irwin, 2015), and ChemSpider databases (Pence and Williams, 2010). Compounds labeled X1, X2, X3, X4, X5, X6, and X7 which showed good half-maximum inhibitory concentration (IC₅₀) values against Leishmania parasite’s sterol methyltransferase together with amphotericin B, miltefosine, paromomycin, and 22,26-azasterol were also used.

The ligands generated from the de-novo design were virtually screened against the LdSMT in the second stage. For both stages, the ligands and the protein were prepared using the AutoDock Tools (Morris et al., 2009) and saved in the input format of AutoDock Vina (Trott and Olson, 2010). The charge, hydrogen bond network, and histidine protonation state of the protein were assigned after pdbqt conversion. Grid box size was set to (91.445 × 73.502 × 78.352) ų with the center at (72.200, 58.009, 13.302) Å. Ligands were then screened against the LdSMT protein with exhaustiveness set to default 8.

The potential protein–ligand complex from the scaffold hopping was submitted to e-LEA3D (Douguet, 2010) for further de-novo design. The binding site radius was set to 15 Å and the final score set to 1 with the same active site coordinates as used in the molecular docking study. Conformational search was set to 10, number of generations to 30, and population size to 30 with the rest of the options left as default.

The atomistic details of binding between the LdSMT and small molecules upon ligand binding were determined using the BIOVIA discovery studio visualizer v19.1.0.18287 (BIOVIA, San Diego, CA, USA) (Šudomová et al., 2019).

The inhibitory constant (K_i ) of the ligands was calculated from the binding energies of the selected compounds and the LdSMT protein (Islam and Pillay, 2020). In addition, ligand efficiency (LE) metrics including ligand efficiency scale (LE_Scale), fit quality (FQ), ligand efficiency-dependent lipophilicity (LELP), and surface binding efficiency were also determined (Hopkins et al., 2014).

The ADMET properties were determined using SwissADME (Daina et al., 2017) and the OSIRIS Property Explorer in Data Warrior (Sander et al., 2015). Pan-assay interference compounds (PAINS) (Daina et al., 2017) and synthetic accessibility (Daina et al., 2017) search using SwissADME (Daina et al., 2017) were performed to eliminate false positive compounds that possess good physiochemical properties as well as those with complex structures.

The biological activity of the selected compounds was predicted using prediction of activity spectra for substance (PASS) (Lagunin et al., 2000) with the simplified molecular input line entry system (SMILES) as inputs.

A 100-ns MD simulation was performed for the unbound LdSMT and the protein–hit complexes using GROMACS 2018 (Van Der Spoel et al., 2005; Abraham et al., 2015). QtGrace (Dahiya et al., 2019) was used to plot the graphs generated from the MD simulation. The binding free energies of the complexes were calculated using MM-PBSA (Kumari et al., 2014). The energy contribution of each residue was also determined using g_MMPBSA. The graphs from the MM-PBSA computations were generated using the R programming language (Tippmann, 2014; Alkarkhi and Alqaraghuli, 2020).

Structural similarity search of all the hits was done via DrugBank (Wishart et al., 2018) to identify drugs with potential antileishmanial activity and possible mechanisms of action from similar compounds.

The 3D structure of LdSMT is yet to be experimentally elucidated; therefore, the structure was modeled. A BLAST (Boratyn et al., 2013) search of the protein sequence of LdSMT (SCMT1/GenBank ID: AAR92099.1) was performed via NCBI BLAST (Boratyn et al., 2013) to identify suitable identical templates to the LdSMT. The search revealed 12 experimentally determined protein structures that are identical to the LdSMT (Supplementary Table 1). The most widely used criteria in selecting a template is to choose the template with the highest sequence identity to the query sequence (Broni et al., 2021) and that was used for the modeling. However, the resolution at which the template protein structure was experimentally determined must also be taken into consideration. Also, the coverage of the template sequence to the query is another important factor. Herein, the E-value, sequence identity, query coverage, and the resolution of the 3D structures were used to select the most suitable templates as previously done (Meier and Söding, 2015; Haddad et al., 2020; Kwofie et al., 2021).

All the 12 identical protein structures had sequence identity less than 30% to the LdSMT, and 5WP4 demonstrated the highest with an identity of 29.01% (Supplementary Table 1); however, 5WP4 had a relatively low coverage of 45% to the LdSMT (Supplementary Table 1) and an E-value of 1 × 10⁻¹¹. The 3BUS template had the least E-value of 6 × 10⁻²⁰ and identity of 24.12%. The 3BUS protein has previously been used in modeling the SMT of L. infantum (Azam et al., 2014). Although 3BUS had a low resolution (2.65 Å), it was selected as one of the structures for modeling. 3BUS is the crystal structure of the rebeccamycin 4′-O-methyltransferase (RebM) in complex with S-adenosyl-l-homocysteine (Singh et al., 2008). On the other hand, the 4PNE template was also shortlisted as a suitable template due to its high coverage to the LdSMT (61%), high resolution (1.50 Å), and sequence identity similar to that of 3BUS (24.15%). 4PNE is the SpnF enzyme in Saccharopolyspora spinosa involved in the biosynthesis of the insecticide spinosyn A (Fage et al., 2015; Jeon et al., 2017). SpnF has been reported to be structurally similar to S-adenosyl-L-methionine (SAM)-dependent methyltransferases (Fage et al., 2015).

Furthermore, a BLAST search via the SWISS-MODEL (Waterhouse et al., 2018) revealed that 4PNE covered residues 43 to 258 while 3BUS spanned from residues 47 to 276 of the LdSMT. From residues 258 to 353 of the LdSMT sequence, both templates do not share similarities with the LdSMT. Thus, the protein structure with PDB ID 6UAK was selected in addition to 3BUS and 4PNE for the multitemplate homology modeling. The 6UAK shared similarity with the LdSMT mostly from residues 100 to 345. 6UAK is a SAM-dependent methyltransferase (LahS_B) from the Lachnospiraceae bacterium C6A11 (Huo et al., 2020). Both 3BUS and 6UAK templates, like the LdSMT, are methyltransferases, while the 4PNE is a methyltransferase-like protein.

An earlier study identified Modeller (Eswar et al., 2008) to predict the most accurate model of L. infantum sterol methyltransferase (Azam et al., 2014). Although the two organisms (L. infantum and L. donovani) belong to the same genus and SMT is highly conserved among Leishmania species (Goto et al., 2007), there was a need to model the structure of LdSMT to ascertain its accuracy. Therefore, Modeller 10.2 (Eswar et al., 2008) was employed for modeling the structure of the LdSMT.

Three modeling approaches were employed to predict the most reasonably accurate LdSMT model. First, 3BUS was used as template to model five structures of the LdSMT. Secondly, 4PNE was also used to model five different structures of the LdSMT. Lastly, a multitemplate homology modeling approach was employed by using three templates comprising 3BUS, 4PNE, and 6UAK protein structures. For each approach, the best model was selected based on the DOPE score, which is an atomic distance-dependent statistical potential calculated from a sample of native protein structures (Shen and Sali, 2006). The DOPE scores were used to distinguish “good” models from “bad” ones with lower DOPE scores signifying a better model (Shen and Sali, 2006; Eswar et al., 2008).

Next, validation and quality assessment of the predicted 3D structures were undertaken to obtain reasonable structures of the proteins. The best models from each of the three different approaches MOD3BUS2, MOD4PNE5, and MOD3TEMP3 were evaluated to select the most reasonably valid structure of the LdSMT.

The percentage of residues in the most favored, additionally allowed, generously allowed, and disallowed regions in a Ramachandran plot determines the quality of protein structures. From the Ramachandran plots generated from PROCHECK, model MOD3BUS2 had 270 (86.8%), 32 (10.3%), 6 (1.9%), and 3 (1.0%) residues in the most favored, additionally allowed, generously allowed, and disallowed regions, respectively (Table 1 and Supplementary Figure 3A). Model MOD3TEMP3 had 266 (85.5%), 31 (10.0%), 8 (2.6%), and 6 (1.9%) residues in the most favored, additionally allowed, generously allowed, and disallowed regions, respectively (Table 1 and Supplementary Figure 3B). For the MOD4PNE5 model, 264 (84.9%) residues were in the most favored region, 32 (10.3%) in the additionally favored region, 11 (3.5%) in the generously allowed region, and 4 (1.3%) in the disallowed region (Table 1 and Figure 3). The Ramachandran plot statistics of all three structures were comparatively close (Table 1) and are consistent with those of a previously modeled LiSMT structure using Modeller (Azam et al., 2014).

The qualities of the three structures were further assessed via SAVES v6.0 (Behera et al., 2021). The protein structure MOD3BUS2 had a VERIFY score of 62.61%, ERRAT quality factor of 41.5663, PROVE score of 11%, and four PROCHECK errors, one warning and three passes (Supplementary Table 3). For the MOD3TEMP3 model, ERRAT predicted an overall score of 45.858, VERIFY a score of 53.82%, PROVE a score of 9.6%, and five PROCHECK errors, one warning and two passes (Supplementary Table 3). Model MOD4PNE5 was predicted to have ERRAT, VERIFY, and PROVE scores of 46.9565, 62.04%, and 7.1%, respectively. MOD4PNE5 was also predicted to have four PROCHECK errors, two warnings and two passes (Supplementary Table 3) For a model to be considered high quality, 80% of its amino acids must have a score of 0.2 in the 3D-1D profile (VERIFY score). Although all the top 3 structures did not have a VERIFY score above 80%, MOD3BUS2 (VERIFY score of 62.61) and MOD4PNE5 (VERIFY score of 62.04) could be considered since a previous study has shown that a crystallized structure also performed poorly based on the VERIFY 3D quality indicator (Mora Lagares et al., 2020). For the ERRAT predictions, MOD4PNE5 had the best score and was also predicted using PROVE to be less erroneous (Supplementary Table 3). Generally, model MOD3TEMP3 exhibited the lowest quality scores from SAVES v6.0.

Based on the quality assessments, the protein structure MOD4PNE5 was selected as the most reasonable structure of the LdSMT protein. Aligning the selected LdSMT structure and the chain A of the 4PNE structure revealed a close similarity with a root mean square deviation (RMSD) of 0.356 Å (Supplementary Figure 2C). The quality of the selected model (MOD4PNE5) was further assessed via ProSA-web (Sippl, 1993; Wiederstein and Sippl, 2007). With a Z-score (Zhang and Skolnick, 1998) of −3.97, the LdSMT structure was predicted to be of X-ray quality (Supplementary Figure 4A). The Z-score was used to indicate the overall protein quality (Zhang and Skolnick, 1998; Benkert et al., 2011). ProSA (Sippl, 1993; Wiederstein and Sippl, 2007) was used to plot the energies as a function of the amino acid sequence position for the LdSMT protein (Supplementary Figure 4B). Positive energy values signify problematic or erroneous regions of the protein structure. The residues of the protein demonstrated relatively low energies until residue 230, where positive energy values were observed until the end of the sequence.

The selected LdSMT structure was then refined using ModRefiner (Xu and Zhang, 2011). The refined structure of the selected LdSMT model showed improved Ramachandran plot statistics with 268 (86.2%), 34 (10.9%), 7 (2.3%), and 2 (0.6%) residues in the most favored, additionally allowed, generously allowed, and disallowed regions, respectively (Table 1).

CASTp was used to predict the plausible binding sites within the refined LdSMT structure. CASTp predicted 65 potential binding sites of the LdSMT protein. Ligand binding sites tend to involve the largest pockets or cavities on the protein (Laskowski et al., 1996; Liang et al., 1998); thus, pockets with relatively low areas and volumes, such that no ligand could fit, were not considered (Broni et al., 2021; Kwofie et al., 2021). A total of seven binding sites were shortlisted (Supplementary Table 4), which were visualized using PyMOL (PyMOL Molecular Graphics System, Version 1.5.0.4, Schrödinger, LLC) (Lighthall et al., 2010) and Chimera 1.16 (Pettersen et al., 2004). However, superimposing the LdSMT on the 4PNE template revealed that pocket 1 was similar to the S-adenosyl-L-homocysteine (SAH) ligand’s binding site in the 4PNE protein (Fage et al., 2015). Surprisingly, pocket 1 for LdSMT was predicted to have no opening in Chimera 1.16 (Pettersen et al., 2004), leaving pockets 2 to 7 as the most plausible binding cavities of the LdSMT (Supplementary Table 4). Pockets 5 and 7 also overlapped and occupied the same region (Supplementary Table 4).

Following the protocols of BRUSELAS, 100 ligands were generated using the ChEMBL database (Davies et al., 2015) with varying degrees of similarities to 22,26-azsterol. Out of the 100 ligands, 17 were identified to possess unique scaffolds devoid of the steroidal nucleus present in the 22,26-azasterol. The total score (comprising a combination of WEGA, LiSiCA, Screen3D, and OptiPharm) ranged from 0.38458 to 0.65814. Ten of the molecules with different scaffolds and varying scores are presented (Supplementary Table 5).

A search via the PubChem (Kim et al., 2021), ZINC15 (Sterling and Irwin, 2015), DrugBank (Wishart et al., 2018), and ChemSpider databases (Pence and Williams, 2010) generated 1,342 derivatives of all 17 scaffolds. The 1,370 compounds consisting of 17 scaffolds, 1,342 derivatives, 22,26-azasterol, 7 other known inhibitors (labeled X1–X7), and 3 already known drugs (amphotericin B, miltefosine, and paromomycin) were screened against an energy minimized LdSMT using a grid box of (91.445 × 73.502 × 78.352) ų with the center at (72.200, 58.009, 13.302) Å to cover the protein. Screening 1,370 compounds against the active site of the protein identified 25 hits which were selected based on the binding affinities and orientation within the binding site of the protein.

Among the three drugs used in leishmaniasis treatment, amphotericin B had the least binding energy of −5.3 kcal/mol followed by paromomycin (−5.0 kcal/mol) and miltefosine (−4.0 kcal/mol). Interestingly, all the known inhibitors had binding energies lower than the three drugs signifying a higher binding affinity to LdSMT. The binding energies of −5.9, −6.2, and −6.5 kcal/mol were obtained for X6, X3, and X4, respectively. The least binding energy of −7.7 kcal/mol was observed for X5 comparable to 22,26-azasterol (−7.6 kcal/mol). In-vitro studies reported that 22,26-azasterol inhibited L. donovani intracellular amastigotes and Trypanosoma brucei subsp. brucei with IC₅₀ values of 8.9 and 1.76 μM, respectively (Magaraci et al., 2003; Gros et al., 2006), supporting the results reported herewith. Compounds X1, X7, and X2 had binding energies of −7.3, −7.2, and −7.0 kcal/mol, respectively (Supplementary Table 6). Similarly, in-vitro studies revealed that X1, X2, X3, X4, X5, X6, and X7 inhibited the growth of Leishmania parasites with IC₅₀ less than 10 μM, except for X6 and X7 which were found to suppress growth with IC₅₀ values of 28.6 and 30 μM, respectively (Magaraci et al., 2003; Lorente et al., 2004; Andrade-Neto et al., 2016; Torres-Santos et al., 2016).

The 12 best hits out of the 25 selected from the scaffold hopping had lower binding energies compared with the three drugs (Supplementary Table 6). In addition, the binding energies of these ligands were also found to be comparable to the known inhibitors with S1 showing the least binding energy of −9.0 kcal/mol. The closest to this ligand were S2, S3, and S4 with binding energies of −8.9, −8.8, and −8.7 kcal/mol, respectively (Supplementary Table 6). Compounds S12, S11, and S10 had the highest binding energies among the 12 best compounds with binding energies of −7.0, −7.0, and −7.2 kcal/mol, respectively. In addition, the binding energies of compounds S9 (−7.3 kcal/mol), S8 (−7.4 kcal/mol), and S7 (−7.4 kcal/mol) were also obtained. Comparable binding energies to the two lowest binding energies of the known inhibitors were obtained for S5 (−7.7 kcal/mol) and S6 (−7.6 kcal/mol).

The de-novo drug design is the generation of novel chemical entities that fit a set of constraints using computational algorithms (Schneider and Schneider, 2016). Despite the challenges of synthetic accessibility associated with this method, the application of de-novo drug design leads to the development of drug candidates in a cost- and time-efficient manner (Mouchlis et al., 2021). In addition, the de-novo design generates novel compounds with improved biological activity (Mouchlis et al., 2021). A number of studies have been undertaken for the de-novo design of inhibitors against plausible targets (Kranthi et al., 2018; Islam and Pillay, 2020; Pathania et al., 2021). A previous study selected a lead molecule based on the least binding energy as well the accurate pose of the ligand within the protein binding pocket for the de-novo design of inhibitors (Pathania et al., 2021). A similar approach was used in the identification of promising anti-DNA gyrase antibacterial compounds (Islam and Pillay, 2020). Ligands with very low binding energies and accurate pose have the potential to inhibit the receptor.

Among the compounds obtained from scaffold hopping, the well-known heterocyclic quinolinone and the phenylpiperazine/phenylpiperidine moieties found in several bioactive compounds were present. These compounds with diverse pharmacological potencies have been explored for various ailments including leishmaniasis (Kshirsagar, 2015; Chanquia et al., 2019; Mishra et al., 2021). One of such derivatives is 2-(4-(4,6-di(piperidin-1-yl)-1,3,5-triazin-2-ylamino)phenyl)-2-methyl2,3-dihydroquinazolin-4(1H)-one which had an IC₅₀ value of 0.65 μM against intracellular amastigotes when compared with miltefosine (IC₅₀ of 8.4 μM) (Sharma et al., 2013), corroborating the fact that de-novo drug design could result in novel scaffolds as potential antileishmanial agents.

The e-LEA3D was used to generate 155 potential novel compounds against LdSMT after using protein–S1, S2, S3, S4, and 22,26-azasterol complexes. They were filtered based on Lipinski’s rule of five (Benet et al., 2016) and redundancy to generate 78 ligands that were subjected to molecular docking studies. The six best hits (A1, A2, A3, A4, A5, and A6) (Figure 4) were selected for downstream analysis. The search via public databases including PubChem (Kim et al., 2021) and ChemSpider (Pence and Williams, 2010) showed that the six ligands do not have duplicates.

The docking analysis of the six compounds revealed the binding energies of A1 (−8.4 kcal/mol), A2 (−7.5 kcal/mol), and A3 (−7.2 kcal/mol), while A4, A5, and A6 had −7.0 kcal/mol. Though A1, A2, A3, and A6 have similar core structures, they had different binding energies (Table 2). Interestingly, using AutoDock Vina, compounds that had shown binding energies ≤−7.0 kcal/mol have been found to demonstrate significant inhibitory activities against the parasite of consideration (Chang et al., 2007; Wyllie et al., 2018; Tabrez et al., 2021). In lieu of this, the predicted compounds may have the potential of suppressing LdSMT since the binding energies were lower than −7.0 kcal/mol. Altogether, compounds A1, A2, A3, A4, A5, and A6 showed binding energies lower than amphotericin B, miltefosine, and paromomycin. Similarly, all the compounds had binding energies comparable to the known inhibitors and, therefore, have the potential of attenuating LdSMT.

Compounds with similar activity against a receptor may possess identical chemical features in sterically consistent locations within the pocket of a macromolecule (Held et al., 2011; Du et al., 2016). Most of the compounds including 22,26-azasterol and its derivatives were observed to interact with residues Asp58, Ala88, Arg89, Tyr92, Glu85, Phe93, Phe100, Glu102, Lys198, Pro199, and Gly200, which lined binding pockets 5 and 7 (Supplementary Tables 4 and 6). A similar study involving L. infantum SMT, however, showed the ligands to interact with residues Tyr1, Gly4, Gln5, Gly45, Gly47, Asn67, Asn68, Gln72, and Ile112 within the binding pocket of the receptor (Azam et al., 2014). The nature of interactions included pi–anion, pi–pi stacking, pi–alkyl, pi–sigma, carbon–hydrogen, and hydrogen bonds similar to the other study (Azam et al., 2014). Among all the ligands, only amphotericin B had five hydrogen bonds, with residues Asp31, Phe307, Val308, Arg309, and Leu310 found to line pocket 2. Compounds X1 and X5 formed two hydrogen bonds each with LdSMT. Compound X1 interacted with pockets 5 and 7 and residues Asp172 and Gly200, while X5 interacted with pocket 4 and residues Asn12 and Thr319 via hydrogen bonds (Supplementary Table 6).

In addition, compounds X2, X6, X7, S1, S5, S9, S10, and S12 docked into the binding pocket of the receptor but showed no hydrogen bond interactions with any of the residues. The compound 22,26-azasterol formed two hydrogen bonds with Glu102 and Gly200 as well as hydrophobic interactions with Phe100, Lys198, and Pro199. Moreover, while paromomycin, S2, S3, S4, S5, and S11 formed a hydrogen bond with Arg89, that of S6 and S7 showed a similar interaction with Asp58. Apart from paromomycin, all the other ligands which interacted with Arg89 had low binding energies, implying that it might be critical for binding. The known drugs and inhibitors formed hydrophobic interactions with one or more of the residues Phe100, Met101, Asp104, Asp172, Pro199, Gly200, Thr201, Tyr343, and Ile344, while the selected S-class compounds showed hydrophobic interactions with one or more of the residues Asp58, Ala88, Arg89, Tyr92, Phe93, and Phe264 (Supplementary Table 6).

For the e-LEA3D-generated hits, apart from A2 and A3 which were predicted not to form hydrogen bond interactions with any amino acid residues, the remaining four exhibited hydrogen bonding with at least one of the amino acid residues of LdSMT (Table 2). This notwithstanding, all six formed hydrophobic interactions with the LdSMT protein. The hydrophobic interactions for A1 were with residues Arg89, Tyr92, Ala95, Ala96, and Leu123 (Figure 5 and Table 2), and those for A2 were with Phe84, Glu85, Ala88, Arg89, and Tyr92 (Table 2 and Supplementary Figure 5A). The ligand A3, on the other hand, formed hydrophobic interactions with Ala88, Arg89, Tyr92, Phe93, Ala96, and Phe264 (Table 2 and Supplementary Figure 5B). Furthermore, both A4 (Supplementary Figure 5C) and A6 (Supplementary Figure 5E) formed hydrogen bonding with Cys202 and Asp58, respectively, while A5 (Supplementary Figure 5D) formed two hydrogen bonding interactions with Arg222 and Lys317.

A druggable candidate must be able to reach the site of action in the body at optimized concentrations (Hughes et al., 2011). Predictions of pharmacological and physicochemical parameters are essential because they offer clues as to whether the molecule could reach the active site in the desired concentration and remain there to elicit the required biological response. The physicochemical and pharmacokinetic parameters of the six compounds were assessed using SwissADME (Hughes et al., 2011). Physicochemical profiling assessed both Lipinski’s rule of five (RO5) and Veber’s rule to determine if the chemical compounds with certain pharmacological or biological activity have chemical and physical properties to make them orally active (Lipinski et al., 2001; Veber et al., 2002). A5 and A6 were found to obey both Lipinski’s and Veber’s rules, and A1, A2, A3, and A4 violated one of the two rules (Supplementary Table 7).

Solubility, an important physicochemical property, was predicted to assess the bioavailability and bioactivity of the hit compounds (van den Anker et al., 2018). Many lead compounds have failed to reach clinical trials despite being potent because of low bioactivity attributed to insufficient solubility, making solubility predictions critical in the early stages of drug design (Das et al., 2022). Compounds A4, A5, and A6 were, however, predicted to be soluble compared with A1, A2, and A3 whose solubility can be improved upon by structural modification (Das et al., 2022).

Next, molar refractivity (MR) was assessed to give valuable information on the pharmacokinetics and pharmacodynamics of the compounds. This is governed by different interactions in solution such as drug–solvent, drug–drug, and drug–co-solute interactions (Sawale et al., 2016). The compounds passed for molar refractivity as the predicted values for MR (Supplementary Table 7) were all within the acceptable range of 40 to 130. Since compounds with topological polar surface area (tPSA) not more than 140 Ų are considered to have good oral bioavailability, all compounds are predicted to be orally active and have good bioavailability.

Furthermore, synthetic accessibility was explored to evaluate the synthetic feasibility of the de-novo hits and has gained importance in the prioritization of compounds in drug design (Ertl and Schuffenhauer, 2009). Many druggable candidates especially de-novo constructed chemical entities are unable to reach clinical trials due to their molecular complexity coupled with difficulty in synthesis (Huang et al., 2010; de Souza Neto et al., 2020). SwissADME (Daina et al., 2017) predicted all the six compounds to possess synthetic accessibility less than 6 implying easy synthesis.

The PAINS of the chemical compounds was investigated to establish whether they will react non-specifically with numerous biological targets (Baell, 2016). A number of PAINS compounds include toxoflavin, isothiazolones, hydroxyphenyl hydrazones, curcumin, phenol-sulfonamides, rhodanines, enones, quinones, and catechols (Baell and Walters, 2014). It is reported that about 5% of US FDA-approved drugs obtained from natural and synthetic drugs still contain PAINS-recognized substructures (Baell and Nissink, 2017). All the compounds were predicted not to contain PAINS substructures.

Pharmacokinetics studies were used to evaluate the time course of absorption, distribution, metabolism, and excretion of the selected hits (Luer and Penzak, 2016). The parameters measured were blood–brain barrier (BBB), gastrointestinal absorption (GI), and permeability glycoprotein (P-gp). Compounds predicted to permeate the BBB have the potential to bind to relevant receptors of the brain to activate signal pathways (Banks, 2009). Four of the compounds (A1, A2, A3, and A5) were predicted not to cross the blood–brain barrier, while A4 and A6 were predicted to cross the BBB to attach to the receptors in the brain to elicit a biological response. GI absorption was probed to investigate whether the hit compounds will be absorbed into the bloodstream after metabolism when orally administered (Löbenberg et al., 2013). All ligands were predicted to have a high GI absorption score except A2 and A3 (Supplementary Table 7). Another pharmacokinetic parameter considered for this study was to explore whether the hits generated were P-gp substrates as compounds predicted to inhibit P-gp result in their increased bioavailability (Lin and Yamazaki, 2003; Prachayasittikul and Prachayasittikul, 2016). The six compounds were screened for their P-glycoprotein binding affinity and they were all predicted to be substrates except A4.

The toxicity profiles of all six compounds were predicted using OSIRIS Property Explorer in Data Warrior (Sander et al., 2015). Toxicity prediction has become very critical in the development of drugs as over 45% of drug candidates fail due to toxicity deficiencies (Van Norman, 2019). Moreover, between 1953 and 2013, as many as 462 medicinal products were withdrawn from the market due to adverse drug reactions (Onakpoya et al., 2016). Toxicity profiling considered for this study was mutagenicity, carcinogenicity, irritancy, and reproductive effects. Of the six compounds, only A1 was predicted to be tumorigenic. The rest neither were mutagenic nor possessed any irritant or reproductive effects (Supplementary Table 8). Among all the compounds under consideration, only X6 was predicted to possess reproductive effects. Overall, the predictions indicate that all molecules may have safe pharmacokinetic and pharmacodynamic profiles except for A1, which would require structural modification to improve its pharmacological properties. For instance, the prediction showed that replacement of the chlorine substituent with a hydrogen atom could render the A1 analog non-tumorigenic.

The open Bayesian machine learning technique, PASS, was used to predict the biological activity of the ligands based on the structure–activity relationship between the selected hits and the training set of compounds of known biological activity (Lagunin et al., 2000; Parasuraman, 2011). A ligand is said to be biologically active and requires experimental validation if the probability of activity (Pa) is greater than the probability of inactivity (Pi) (Basanagouda et al., 2011). Among the six compounds, A3 was predicted to possess antileishmanial properties with a Pa of 0.362 and a Pi of 0.066 and also dermatological properties with a Pa of 0.32 and a Pi of 0.091. Compounds A5 and A6 were also predicted as dermatologic, with Pa values of 0.205 and 0.249 and Pi values of 0.162 and 0.120, respectively. The results may suggest that A3, A5, and A6 might be beneficial in treating post-kala-azar leishmaniasis (Momeni et al., 2003; Ali et al., 2012).

Compounds A1 (Pa of 0.571 and Pi of 0.111), A3 (Pa of 0.774 and Pi of 0.026), and A6 (Pa of 0.615 and Pi of 0.079) were also predicted as mucomembrane protectors. A recent in-vitro study revealed that butein acting as a mucomembrane protector on human cells increased immunity against pathogenic infections (Satari et al., 2021). This may suggest that the compounds have the potential of boosting the immune system to prevent disease exacerbation. Compound A1 was predicted as an indolepyruvate C-methyltransferase inhibitor with a Pa of 0.226 and a Pi of 0.074. Compound A3 was also predicted to be phenol O-methyltransferase, histamine N-methyltransferase, and acetylserotonin O-methyltransferase inhibitors with Pa values greater than Pi.

The inhibitory constant (K_i ) and other parameters such as LE, LE_Scale, FQ, binding efficiency index (BEI), surface efficiency index (SEI), and LELP were calculated (Supplementary Table 9).

K_i is the concentration required to produce half-maximum inhibition and, hence, an indicator of the potency of a ligand (Fisar et al., 2010). Computation of K_i for the protein–ligand complexes was obtained using Equation (1), where R is the molar gas constant (1.987 × 10⁻³ kcal/K mol⁻¹) and T (298.15 K) is the absolute temperature (Du et al., 2016).

The K_i predicted for the ligands was low (Supplementary Table 9), hence has the capacity to be lead-like with possible high potency (Reynolds and Reynolds, 2017).

LE is a value that expresses the binding energy of a compound normalized by the compound’s size and expressed by the number of heavy (non-hydrogen) atoms (Hopkins et al., 2004). This property is important to consider in screening for hits as larger compounds tend to show greater binding energy due to a large number of interactions but may not necessarily be the most efficient binder (Hevener et al., 2018). The LE was computed using Equation (2), where BE is the binding energy and NHA is the number of heavy atoms (Abad-Zapatero et al., 2010).

Interestingly, all the ligands except A1, A3, and A5 were within the optimal range (LE < 0.3 kcal/mol/HA) (Schultes et al., 2010) for the ligand efficiency of lead-like molecules.

Results from the computation of LE being size-dependent may not be a true reflection of the binding energy of the compound, and therefore, ligand efficiency scaling (LE_Scale), a size-independent parameter that compares ligands with the help of an exponential function to the maximal LE values, is required (Reynolds et al., 2007). LE_Scaling was computed using Equation (3), where NHA is the number of heavy atoms (Reynolds et al., 2007).

Potential lead-like molecules are suggested to have an LE_Scale lower than 0.3 (Islam and Pillay, 2020). The LE_Scale values of all six molecules were, however, predicted to be in the range of 0.3 to 0.5 with A5 having the highest LE_Scale of 0.455.

FQ is another size-independent parameter that determines the optimal ligand binding within the receptor active site (Schultes et al., 2010), and is computed using Equation (4) (Schultes et al., 2010).

FQ scores range from 0 to 1 with values close to 1 signifying an optimal ligand binding (Schultes et al., 2010). All the compounds were predicted to have an FQ score above 0.7 except A6 (0.695), implying a stronger ligand binding. With the predicted FQ being close to 1, it suggests an optimal ligand binding.

The LELP, on the other hand, assesses the binding energy of a compound in relation to the compound’s lipophilicity (Hevener et al., 2018). LELP is a parameter used in drug design and development to evaluate the quality of compounds by linking potency and lipophilicity in an attempt to estimate drug-likeness (Edwards and Price, 2010). Equation (5), where logP is the lipophilicity, was used in calculating the LELP of the compounds (Schultes et al., 2010).

The recommended range for promising molecules for LELP was >3 (Schultes et al., 2010). All the analogs had LELP above 4 (Supplementary Table 9) suggesting an optimized affinity with respect to lipophilicity.

BEI and SEI are two alternative metrics that are also used to compare the activity of molecules according to size and area (Schultes et al., 2010). Binding efficiency index is defined by BEI = p(IC₅₀)/MW, where MW is the molecular weight (Schultes et al., 2010). The relation in Equation (6) was used in computing the BEI of the ligands (Abad-Zapatero et al., 2010).

SEI, on the other hand, is defined by SEI = p(IC₅₀)/PSA, where PSA is the polar surface area of the ligand (Abad-Zapatero et al., 2010). Calculation of SEI was done using Equation (7) (Abad-Zapatero et al., 2010).

By rule of thumb, potential inhibitors must approximately have the same BEI and SEI values (Abad-Zapatero, 2007). Compounds A1, A4, and A6 were predicted to have BEI equal to SEI. Altogether, the parameters predicted for all the compounds were mostly within the acceptable range prompting the need for experimental analysis.

MD simulation is a computer simulation method for analyzing the physical movements of atoms and molecules (Hollingsworth and Dror, 2018). Of great concern for the MD simulation is how a biomolecular system responds to some perturbation within a short period of time (Karplus and Kuriyan, 2005; Hollingsworth and Dror, 2018). Due to its usefulness, a number of drug design studies have explored molecular dynamics simulation to analyze and validate the binding poses, stability of the complexes, and binding affinity of selected hits within the binding pocket of the receptor (Agyapong et al., 2021; Gupta et al., 2021; Pathania et al., 2021). To check the relative stability of each complex using a 100-ns time span MD simulation, parameters such as RMSD, root mean square fluctuation (RMSF), and radius of gyration (Rg) were computed.

Four (A1, A2, A3, and A6) out of the six compounds possessed benzo[b]azepine moiety as a replacement for the steroidal core in the 22,26-azasterol. This same moiety is present in paullone and its derivatives as well as BNZ-1 (Figure 8), which are known to suppress growth in Leishmania parasites with IC₅₀ values of 47 nM and 100 μM, respectively (Clark et al., 2007; Dao Duong Thi et al., 2009), implying the possible antileishmanial potentials of the proposed hits. Moreover, the chemical structural similarity search of compounds A1, A2, A3, and A6 via DrugBank (Wishart et al., 2018) revealed a variable similarity to antipsychotic and antidepressant drugs (Kaur, 2013; Mendonca Junior et al., 2015). For instance, compounds A1, A2, A3, and A6 showed similarity scores above 0.50 to vabicaserin, sertraline, and indatraline (Richardson et al., 2009; Wishart et al., 2018; Lima et al., 2018). In addition, the four compounds showed similarity scores of around 0.55, which are close to those of daledalin, zanapezil, clocapramine, imipramine, and dimethacrine (Figure 8) (Mukherjee et al., 2012; da Silva Rodrigues et al., 2019). Interestingly, these drugs have been explored for their antileishmanial potentials causing Leishmania parasites to undergo mitochondrion depolarization in addition to inhibiting trypanothione reductase, thereby inducing strong oxidative stress in the parasite (Kaur, 2013; Andrade-Neto et al., 2016; da Silva Rodrigues et al., 2019). The similarity scores and the antileishmanial properties of these antidepressant and antipsychotic drugs warrant the testing of the compounds to assess their antileishmanial propensity.

Furthermore, compound A4 devoid of the benzo[H]azepine also showed a chemical structural similarity score of 0.564 to CP-39,332, a serotonin–norepinephrine reuptake inhibitor. The successes emanating from other serotonin inhibitors for leishmaniasis treatment suggest A4 as a potential antileishmanial compound. In addition, A5 showed a similarity score of 0.538 to cefradine, a broad-spectrum antibiotic for the treatment of skin, chest, throat, and ear infections (Wishart et al., 2018). A recent study has revealed that patients exposed to antibiotics had a greater healing rate (Barakat et al., 2017), suggesting A5 to be explored as an antileishmanial agent. In lieu of the aforementioned, the potential leads A1, A2, A3, A4, A5, and A6 with diverse structural similarities with the antipsychotic and antibiotic agents are also proposed as potential antileishmanial agents via inhibition of sterol methyltransferase and are worthy of further experimental evaluation to assess their biological efficacy.