Computational study to discover potent phytochemical inhibitors against drug target, squalene synthase from Leishmania donovani.

AIMS: The parasite, Leishmania donovani is responsible for lethal visceral leishmaniasis (VL) in humans. There is a need to investigate novel medicines as antileishmanial drugs, as medication currently introduced for leishmaniasis may cause resistance, serious side-effects, chemical instability and high cost. Therefore, this computational study was designed to explore potential phytochemical inhibitors against Leishmania donovani squalene synthase (LdSQS) enzyme, a drug target. MAIN
METHODS: Multiple sequence alignment was carried to detect conserved regions across squalene synthases from different Leishmania spp. Their evolutionary relationships were studied by generating phylogenetic tree. Homology modeling method was used to build a three dimensional model of the protein. The validated model was explored by docking simulation with the phytochemicals of interest to identify the most potent inhibitors. Two reported inhibitors were used as references in the virtual screening. The top hit compounds (binding energy less than -9 kcal/mol) were further subjected to intermolecular interaction analysis, pharmacophore modeling, pharmacokinetic and toxicity prediction. KEY
FINDINGS: Seven phytochemicals displayed binding energies less than -9 kcal/mol hence demonstrating ability to be strongly bound to the active site of LdSQS to inhibit the enzymatic activity. Ancistrotanzanine B demonstrated the lowest binding affinity of -9.83 kcal/mol superior to reported inhibitors in literature. Conserved two aspartate rich regions and two signatory motifs were found in the L. donovani squalene synthase by multiple sequence alignment. In addition, study of pharmacophore modeling confirmed that top hit phytochemicals and the reported inhibitor (E5700) share common chemical features for their biochemical interaction with LdSQS. Among seven phytochemicals, 3-O-methyldiplacol showed admissible physicochemical, pharmacokinetic and toxicity predictions compared to the reported inhibitors. All seven phytochemicals satisfied in silico prediction criteria for oral bioavailability. SIGNIFICANCE: Based on the current study, these hits can be further structurally optimized and validated under laboratory conditions to develop antileishmanial drugs.

Introduction

The parasite, L. donovani is one of the causative agents for deadly visceral leishmaniasis (VL) which is transmitted to humans by the bite of infected sand flies. It severely damages the human internal organs becoming highly lethal [1].

The front-line chemotherapeutics [2, 3, 4, 5] for VL are not very efficient to eradicate the disease majorly due to high cost, drug resistance and clinical side effects [6]. Hence it is urgent to investigate new medications that are inexpensive and safe to treat VL.

The L. donovani squalene synthase (LdSQS) is the first enzyme participating in the reaction chain of sterol metabolism in the parasite. The enzyme catalyzes the dimerization of two farnesyl pyrophosphate molecules to produce squalene which then undergoes different enzymatic reaction steps to become ergosterol in Leishmania parasites in contrast to cholesterol in humans [7]. Ergosterol is the abundant membrane sterol in Leishmania and Trypanosoma parasites, which is essential for their growth and proliferation [8]. There have been studies of potent inhibitors of squalene synthase [9, 10], squalene epoxidase [11], C-14 alpha demethylase [12, 13], and sterol-24 methyltransferase [14, 15] in the sterol synthesis pathway, with varying degrees of success.

Zaragozic acid A and quinuclidine derivatives are known inhibitors against squalene synthase from Leishmania spp and Trypanosoma parasites. Two derivatives of quinuclidine, ER119884 and E5700 have shown efficient inhibition at subnanomolar concentrations under long incubation periods against L. amazonensis squalene synthase [10].

The inhibition of LdSQS would lead to lysis of the parasite due to the depletion of endogenous sterol levels, thereby becoming a chemotherapeutic target to investigate potent inhibitors. The LdSQS shared only a 39% identity with the human squalene synthase after the pairwise sequence alignment using EMBOSS Needle server [16] (https://www.ebi.ac.uk/Tools/psa/emboss_needle/). This reveals that the parasite target is substantially different from the human homolog thus can be a promising drug target to discover antileishmanial drugs.

A pilot docking study using homology modeled LdSQS with the reported inhibitor (E5700) and two plant derived alkaloids (ancistrotanzanine b and holamine) was conducted and published by the first author [17]. The promising results from the pilot study led to design a full scale computational study to explore potent inhibitors from plant derived compounds against LdSQS.

In the current study, the modeled LdSQS protein was subjected to further energy minimization using steepest descent algorithm than in the previously published model [17]. This brought significant improvement in the ERRAT and MolProbity validation results of the modeled protein thus enhancing the accuracy of the new LdSQS model. Then the constructed compound library was used to discover potent squalene synthase inhibitors by molecular docking using AutoDock tool [18]. The top hit compounds were further evaluated using pharmacophore modeling and in silico pharmacokinetic and toxicity predictions. To the best of our knowledge, this is the first full scale computational study of modeled LdSQS to identify potent phytochemical inhibitors.

Materials and methods

Multiple sequence alignment, conserved domain analysis, motif identification and phylogenetic tree construction

The squalene synthase sequences (FASTA format) of seven Leishmania species were retrieved from the Universal Protein Knowledgebase (UniProtKB) database (https://www.uniprot.org/) [19]. A multiple sequence alignment for these protein sequences was carried out using Clustal Omega server (https://www.ebi.ac.uk/Tools/msa/clustalo/) [20] and represented by ESPript 3.0 server [21]. The conserved domain database [22] was utilized to identify conserved domains present in the squalene synthase sequences whereas the MOTIF search server (https://www.genome.jp/tools/motif/) was exploited to find conserved signatory motifs using PROSITE database. The phylogenetic tree was constructed using One Click option by Phylogeny. fr server [23].

Homology modeling and validation

Based on the methodology followed in the previously published pilot study, 50 models of the LdSQS protein was built [17]. The solved crystal structure of LdSQS is not available at the RCSB Protein Data Bank (https://www.rcsb.org/) [24]. The amino acid sequence of LdSQS was retrieved from the UniProtKB database in FASTA format (accession no: Q257D4). The BLASTp tool was used to select the best homologous protein template based on sequence identity, query coverage and E value. The final pairwise sequence alignment of template and target was carried out in Clustal Omega server. Prior to the final pairwise sequence alignment, terminal segments of both ends in the target sequence were eliminated as these residues are absent in the template. Modeller v9.22 [25] was used to build 50 models of LdSQS. The model with the least discrete optimized protein energy (DOPE) score was selected for further evaluation. Then the selected model was refined by quality information from MolProbity server [26] using UCSF Chimera v.1.12 [27]. The selected model was energy minimized up to 2100 steps by steepest descent algorithm in AMBER force field (using UCSF Chimera) instead of 1000 steps which were carried in the previous publication [17]. This additional step not only relieved steric clashes and improper contacts of the model but also that increased the values of validation results hence improving the quality of the model. Eventually the quality of optimized model was further assessed using MolProbity, SAVES v5.0 package (https://servicesn.mbi.ucla.edu/SAVES/), ProSA [28] and QMEAN [29] web platforms. The most accurate model of LdSQS was employed for active site prediction and virtual screening with the compounds library. Physicochemical properties and protein secondary structure analysis of the modeled protein were performed by ProtParam [30] and PDBsum [31] web servers respectively.

Identification of binding site residues

The prominent binding site residues of LdSQS were predicted based on the following criteria.

Referring the corresponding publication of selected template (PDB ID: 3WCA) [9] used for homology model building

Selecting residues within 5 A° of the reported inhibitor 1 (E5700) (taken from PDB ID: 3WCC) after superimposition of LdSQS modeled protein with co-crystalized protein template (PDB ID: 3WCC)

Docking of reported inhibitor 1 (E5700) with LdSQS modeled protein

Construction of compounds library

Phytochemicals have diverse and complex structures of scaffolds which can be used as chemical templates for rational drug design. An antiparasitic phytochemical library was prepared using 31 compounds comprising naphthylisoquinoline alkaloids [32, 33, 34, 35, 36, 37, 38, 39, 40], monoterpenoid indole alkaloids [41], c-geranyl flavonoids [42], (+)-trans-hexahydrodibenzopyrans [43], benzoquinolizidine alkaloids [44], anthranoid compounds [45], acridone alkaloid derivatives [46], coumarins [47], benzophenones [48], aporphine alkaloids [49, 50]and pyrimidine-ß-carboline alkaloids [51] to be explored for their potential antileishmanial activity. Two reported squalene synthase inhibitors (E5700 & ER119884) were used as references throughout the computational study.

Grid box preparation and virtual screening

Grid box parameters were prepared covering all the binding site residues and giving optimum sized cavity space to accommodate each ligand during molecular docking process.

The hydrogens were added to ligands and homology modeled protein (LdSQS) at physiological pH (7.4) using Avogadro v1.2.0 [52] software. All phytochemicals were then subjected to addition of gastieger charges followed by energy minimization prior to molecular docking. All histidine residues of LdSQS were kept as neutral during the molecular docking process. Molecular docking calculations were performed by AutoDock Tools 1.5.6. Rigid docking simulation was carried out in Lamarckian Genetic Algorithm (LGA) with a population size of 150 with 2 500 000 evaluations and 27 000 generations. Default values were used for all other parameters. The modeled LdSQS was kept rigid while ligands were flexible throughout the rigid docking method. The ligand efficiency was calculated by, LE = -ΔG/HA, where -ΔG denotes the binding affinity of the ligand and HA is the number of non-hydrogen atoms present in the ligand [53]. In addition, the three dimensional (3D) and two dimensional (2D) interaction diagrams of the best docking complexes were visualized using UCSF Chimera and BIOVIA Discovery Studio Visualizer 2020 Client [54] respectively.

Pharmacophore modeling

The best ligands (compounds) that exhibited binding affinities less than -9 kcal/mol (threshold value), were analyzed for pharmacophoric features using PharmaGist web server [55]. Thereafter, the output files obtained by PharmaGist were submitted to ZINCPharmer server [56] to visualize the best pairwise structural alignment of ligands (phytochemicals) with the reported inhibitor, E5700 (pivot molecule). Scores were calculated using default feature weighting values assigned to each pharmacophore feature by PharmaGist server.

Physicochemical properties, drug-likeness and in silico pharmacokinetic and toxicity prediction

The physicochemical properties of selected ligands were studied using the Molinspiration server (https://www.molinspiration.com/) [57]. Lipinski's rule of five (LRo5) [58] and Veber rule [59]were applied to predict drug-likeness of selected ligands to explore their oral bioavailability. In the present study, absorption, distribution, metabolism, excretion and toxicity (ADMET) properties of the compounds were predicted using pkCSM server [60].

The pharmacodynamics of miltefosine and top hit phytochemicals were explored using available literature based on the half maximal inhibitory concentration (IC50) against Leishmania spp. and Trypanosoma parasite. Further, the pharmacodynamics of E5700 and ER119884 were investigated through literature search in terms of IC50 values against squalene synthase of Leishmania amazonensis.

Results and discussion

Multiple sequence alignment revealed most of the amino acid residues are conserved throughout the squalene synthase sequences (SQSs) of Leishmania species. The conserved domain analysis confirmed SQSs belong to the same isoprenoid biosynthesis enzymes, class 1, superfamily with substrate binding pockets, substrate-Mg2+ binding sites, active site lid residues, catalytic residues and aspartate-rich regions 1 and 2. All amino acid residues belonged to those domains were highly conserved across seven SQSs. The multiple sequence alignment is illustrated in Figure 1B. The aspartate rich regions in LdSQS were positioned at 72DTLED76 and 212DYYED216. Similarly, previous studies reported aspartate rich motifs in human SQS (80DTLED84 & 219DYLED223) [61], Panax ginseng SQS (77DTVED81 & 213DYLED217) [62], cucurbitaceae SQS (77DTVEDD82 & 213DYLED217) [63] and solanaceae SQS (40DTVED44 & 180DYLED184) [64]. These aspartate rich motifs are binding sites for Mg2+ ions. The co-factor (Mg2+) plays a key role in fixing and activating the diphosphate moiety of substrate to form an allylic carbocation during condensation reaction [65]. Our homology modeled protein (LdSQS); Mg2+ interacts with Asp72, Asp76 and Glu75 residues of 72DTLED76 region.

Figure 1

A) Phylogeny of squalene synthase sequences from seven Leishmania species, L. guyanensis squalene synthase (LgSQS), L. panamensis squalene synthase (LpSQS), L. braziliensis squalene synthase (LbSQS), L. infantum squalene synthase (LiSQS), L. major squalene synthase (LmSQS), L. donovani squalene synthase (LdSQS) and L. mexicana squalene synthase (LmeSQS). Two separate clusters are enclosed in blue colour boxes. The phylogenetic tree was constructed by Phylogeny. fr.server. B) Multiple sequence alignment of squalene synthase sequences from Leishmania species. Identical and similar amino acid residues are highlighted with red and yellow background colors. The conserved two aspartate rich regions are presented in green colour amino acid, one letter codes. The purple color triangles and blue color circles under the sequences indicate the conserved catalytic site residues and substrate binding pocket residues respectively. The conserved signatory motif 1 and 2 are highlighted in blue and black color boxes. The figure was generated by ESPript 3.0 tool. The graphical representation of amino acid sequence conservation was created using the WebLogo tool.

Motif Search tool identified two conserved motifs (16 residues long, 164YCHYVAGLVGHGLTQL179 and 28 residues long, 200MGLFLQKTNIIRDYYEDIREEPPRMFWP227) with the conserved signatory motif patterns of Y-[CSAM]-x(2)-[VSG]-A-[GSA]-[LIVAT]-[IV]-G-x(2)-[LMSC]-x(2)-[LIV] and [LIVM]-G-x(3)-Q-x(2,3)-[ND]-[IFL]-x-[RE]-D-[LIVMFY]-x(2)-[DE]-x(4,7)-R-x-[FY]-x-P. These conserved two motif sequences were previously reported in SQS derived from cucurbitaceae [63] and solanaceae [64] family plants with few residual substitutions.

The phylogenetic tree (Figure 1A) was constructed to understand the evolutionary relationships between LdSQS and other SQSs from different Leishmania spp. deposited in UniProtKB Database. Mainly two clades were observed in the phylogeny. One clade consisted of four squalene synthases while other contained three. The phylogenetic tree showed LdSQS is more closely related to L. mexicana SQS followed by L. major and L. infantum. The other clade consisted of three sequences showed the most distant relationship regarding LdSQS.

The protein sequence of LdSQS (comprised of 414 amino acids) was retrieved from UniProtKB Database. The best homologous template (PDB ID: 3WCA_A solved at 2.24 A° resolution) showed 81% query coverage, 61.06% sequence similarity and a 5e−153 E-value with the target protein sequence confirming a single template can be used for model building. A low E-value always indicates high protein sequence similarity between query sequence and template. Generally, a single template can be used to build a quality homology model of target protein when its sequence similarity is more than 60%.The Clustal Omega tool was used to perform the pairwise sequence alignment of target protein and template. The alignment indicates the presence of conserved aspartate-rich regions (43DTLED & 183DYYED187) and binding site residues. However both C and N termini segments of LdSQS sequence were removed before the optimal alignment (Figure 2C). These removed segments have no role on binding site and catalytic activity of the LdSQS protein.

Figure 2

A) Structure of homology modeled protein, LdSQS. The N-terminal and C-terminal are represented in blue colour and red colour respectively. B) The superimposition of template protein 3WCA_A (Yellow) and the 3D modeled structure of LdSQS (Green) (RMSD = 0.308A°). C) The pairwise sequence alignment of LdSQS (Target) with the template protein PDB ID: 3WCA_A. Identical and similar amino acids are highlighted with red and yellow background colors. Dissimilar amino acids are represented in white background. The conserved aspartate rich regions are highlighted in green colour, one letter codes. The symbols, η refers to a 310-helix. The α helices, 310-helices and π helices are displayed as medium, small and large squiggles respectively. Strict ß turns are shown as TT letters. This figure was prepared using ESPript 3.0. The logos of pairwise sequence alignment were generated by WebLogo tool.

Initially 50 structural models of LdSQS were built using the template structure of putative, farnesyltransferase of Trypanosoma cruzi strain (PDB ID: 3WCA_A) (Table S1A). The well optimized homology model was achieved by following the procedure mentioned in the materials and methods section. Eventually the top modeled protein was subjected to quality validation using MolProbity [26], PROCHECK [66], ProSA [28], ERRAT [67], Verify3D [68] and QMEAN [29] web platforms.

The root mean square deviation (RMSD) was calculated by superimposing homology model and template on each other. The RMSD value was found to be 0.308 A° (Figure 2B) suggesting close relationship between the two structures. Details of structure validation results are summarized in Table S1B and illustrated in Figure S1.

The LdSQS is a monomeric protein which predominantly consists of alpha helices and beta turns. The predicted model was then deposited in the protein model database with PMDB ID: PM0083465 [69]. The 3D structure of LdSQS protein is depicted in Figure 2A. The predictions of physicochemical properties and protein secondary structure analysis of the LdSQS protein are shown in Table 1 and Figure 3 respectively.

Table 1

Physicochemical properties of LdSQS (homology model) by ProtParam tool.

ProtParam Parameters	Values
Number of amino acids	414
Molecular weight	47344.59 kDa
Theoretical pI	5.65
Percentage amino acid composition	Ala-28 (6.8%), Arg-25 (6.0%), Asn-12 (2.9%), Asp-31 (7.5%),
	Cys-14 (3.4%), Gln-17 (4.1%), Glu-25 (6.0%), Gly-17 (4.1%),
	His-7 (1.7%), Ile-20 (4.8%), Leu-39 (9.4%), Lys-24 (5.8%),
	Met-17 (4.1%), Phe-20 (4.8%), Pro-16 (3.9%), Ser-29 (7.0%),
	Thr-18 (4.3%), Trp-4 (1.0%), Tyr-17 (4.1%), Val-34 (8.2%)
Atomic composition	Carbon (C) = 2103, Hydrogen (H) = 3302, Nitrogen (N) = 560,
	Oxygen (O) = 620, Sulfur (S) = 31
Formula	C2103H3302N560O620S31
Total number of atoms	6616
Number of negatively charged residues	(Asp + Glu) = 56
Number of positively charged residues	(Arg + Lys) = 49
Extinction coefficient	48205 [Abs 0.1% (=1 g/l)] 1.018, assuming all pairs of Cys residues from cysteines47330 [(Abs 0.1% (=1 g/l) 1.000, assuming all Cys residues are reduced
Estimated half-life	30 h (mammalian reticulocytes, in vitro)
	>20 h (yeast, in vivo)
	>10 h (Escherichia coli, in vivo)
Instability index	38.88
Aliphatic index	86.16
Grand average of hydropathicity (GRAVY)	-0.157

Figure 3

Protein secondary structure analysis of the LdSQS protein using PDBsum tool.

The important binding site residues were identified based on the three step methodology mentioned in the materials and methods section. All identified binding site residues are given in Table S2.

The compounds library was constructed using thirty one phytochemicals. These compounds have previously demonstrated in vitro antiparasitic activities against Leishmania and Trypanosoma parasites. In total, 28 compounds have been reported in vitro antileishmanial activity excluding ancistrobertsonine A, ancistrotanzanine C and ancistrocongoline A. But their exact mechanism of action is undisclosed. Structures of phytochemicals were obtained by the PubChem CID into UCSF Chimera and saved as pdb files except for 5-epi-4′-O-demethylancistrobertsonine C. The structure of 5-epi-4′-O-demethylancistrobertsonine C was first drawn using ACD/ChemSketch Freeware version and then converted to simplified molecular-input line-entry system (SMILE) notation. Then the SMILE notation was used to generate the pdb file using the structure building option in UCSF Chimera. Both reported inhibitors (E5700 and ER11884) were taken from their co-crystalized protein structures, PDB ID: 3WCC and PDB ID: 3WCB respectively.

A library containing 33 compounds was constructed for virtual screening with the homology modeled protein. All the 2D structures of 33 compounds are depicted in Figure 4.

Figure 4

2D chemical structures of the phytochemicals and reported inhibitors used for this study.

AutoDock 4.2 was utilized to perform molecular docking which can be used to identify the binding affinities and specific intermolecular interactions between ligands and protein. The grid box parameters selected for molecular docking study are given in Table S2. All the AutoDock calculations are tabulated in Table 2. The docked ligands (compounds) exhibited binding affinities in the range of -6.87 kcal/mol to -9.83 kcal/mol and their inhibition constants were between 61.84 nM-9.28 μM range. The inhibition constant (Ki) is inversely proportional to the binding affinity of the ligand to the enzyme. Higher the value of binding affinity lower is the value of inhibition constant. The ligands with higher binding affinities are strongly bound to the active site of the protein. Experimentally known two inhibitors were used to compare docking results. Both inhibitors are quinuclidine derivatives namely E5700 (reported inhibitor 1) and ER119884 (reported inhibitor 2). Out of 31 phytochemicals, only one compound showed higher docking score than E5700 whereas 25 compounds exhibited higher docking scores than ER119884. Consequently, -9.00 kcal/mol was set as the threshold level of binding affinity.

Table 2

Docking results of phytochemicals and reported inhibitors used for this study.

No	Ligand Name	PubChem CID	LBE∗ (kcal/mol)	EIC∗∗	LE∗∗∗	Plant Source	Ref.
1	Ancistrotanzanine B	10002166	-9.83	61.84nM	0.32	A	[32]
2	E5700	56947067	-9.81	64.82nM	0.31	N/A#	[9]
3	Ancistrotanzanine A	10319100	-9.68	80.37nM	0.32	A	[32]
4	Ancistrobenomine A	11464259	-9.54	101.91nM	0.31	B	[33]
5	5-epi-4′-O-demethylan cistrobertsonine C	SMILE##	-9.34	143.41nM	0.30	C	[34]
6	Ancistroealaine B	10364043	-9.32	146.86nM	0.31	D	[35]
7	Ramiflorine B	101691176	-9.23	171.07nM	0.26	E	[41]
8	3-O-methyldiplacol	21607150	-9.10	213.00nM	0.28	F	[42]
9	Ancistrogriffine C	11003736	-8.97	263.97nM	0.31	G	[36]
10	Ancistrobertsonine A	15840195	-8.97	264.57nM	0.29	H	[37]
11	Ancistrotanzanine C	11350305	-8.91	293.00nM	0.30	A	[38]
12	Ancistrocongoline C	636651	-8.90	301.44nM	0.29	C	[39]
13	Machaeridiol C	10882982	-8.85	328.04nM	0.33	I	[43]
14	Ancistrocongoline B	11742933	-8.75	386.51nM	0.28	C	[39]
15	4′-O-methyldiplacone	24854122	-8.73	397.43nM	0.27	F	[42]
16	3′-O-methyldiplacone	14539951	-8.68	430.55nM	0.27	F	[42]
17	Ancistrocladidine	193880	-8.66	450.72nM	0.29	J	[40]
18	Cephaeline	442195	-8.61	492.04nM	0.25	K	[44]
19	Ancistrogriffine A	101204923	-8.53	555.63nM	0.29	G	[36]
20	Ancistrocongoline A	392424	-8.49	599.03nM	0.29	C	[39]
21	Vismione D	5281573	-8.40	692.65nM	0.28	L	[45]
22	Ancistrocongoline D	10409129	-8.39	703.8nM	0.28	C	[39]
23	Korupensamine A	392421	-8.35	759.86nM	0.30	C	[39]
24	Ancistroheynine B	135508453	-8.13	1.09μM	0.28	J	[40]
25	Gravacridonediol	5317836	-8.03	1.31μM	0.32	M	[46]
26	Machaeridiol B	10361425	-8.03	1.29μM	0.30	I	[43]
27	ER119884	9844255	-8.02	1.32μM	0.27	N/A#	[9]
28	Mammea A/BB	11750116	-7.79	1.95μM	0.26	N	[47]
29	3-Geranyl-2, 4, 6-trihydroxybenzophenone	6477676	-7.74	2.13μM	0.29	O	[48]
30	Klugine	10457340	-7.64	2.52μM	0.22	K	[44]
31	Cryptodorine	11438278	-7.34	4.13μM	0.32	P	[49]
32	Annomontine	5257090	-7.22	5.11μM	0.36	Q	[51]
33	Duguetine	45102748	-6.87	9.28μM	0.26	R	[50]

LBE∗- Ligand Binding Energy, EIC∗∗- Estimated Inhibition Constant, LE∗∗∗- Ligand Efficiency.

N/A#- Not Applicable.

SMILE##- SMILE notation (Cc1cc2c(c(O)c1)c(OC)ccc2c1c(cc(OC)c2[C@@H](C)N(C)[C@@H](C)Cc12)OC)

A- Ancistrocladus tanzaniensis, B- Ancistrocladus benomensis, C- Ancistrocladus congolensis. D- Ancistrocladus ealaensis, E- Aspidosperma ramiflorum, F- Mimulus bigelovii, G- Ancistrocladus griffithii, H- Ancistrocladus robertsoniorum, I- Machaerium multiflorum, J- Ancistrocladus heyneanus, K- Psychotria klugii, L- Vismia orientalis, M- Thamnosma rhodesica, N- Calophyllum brasilense, O- Garcinia vieillardii, P- Guatteria dumetorum, Q- Annona foetida, R- Duguetia furfuracea.

Eight ligands were identified surpassing the threshold value. Hence they were used for further investigation. In the top hit ligands, ancistrotanzanine B showed the best docking pose with the least binding energy (-9.83 kcal/mol and Ki of 61.84 nM) with LdSQS closely followed by reported inhibitor 1 (E5700), ancistrotanzanine A, ancistrobenomine A, 5-epi-4′-O-demethylancistrobertsonine C, ancistroealaine B, ramiflorine B and 3-O-methyldiplacol with binding affinities of -9.81 kcal/mol (Ki of 64.82 nM), -9.68 kcal/mol (Ki of 80.37 nM), -9.54 kcal/mol (Ki of 101.91 nM), -9.34 kcal/mol (Ki of 143.41 nM), -9.32 kcal/mol (Ki of 146.86 nM), -9.23 kcal/mol (Ki of 171.07 nM) and -9.10 kcal/mol (Ki of 213.00 nM) respectively. Interestingly, the docking scores for E5700 and ancistrotanzanine B in the current study, were slightly different than the previously published results where E5700 (-9.75 kcal/mol) exhibited higher score than ancistrotanzanine B (-9.55 kcal/mol) [17].

Among the top hit ligands, except ramiflorine B and 3-O-methyldiplacol, the other ligands exhibited ligand efficiency greater than 0.29 kcal/mol/HA which is considered as good starting hits for lead optimization. Instead of considering the affinity of the whole compound the average affinity per atom is examined during ligand efficiency calculations.

Seven phytochemicals and reported inhibitor 1 (E5700) were recognized for further analysis as they displayed strong docking scores than the threshold level.

Interestingly, ancistrotanzanine B, ancistrotanzanine A and ramiflorine B showed no conventional hydrogen bonds. Instead they were stabilized via hydrophobic interactions at the binding site of LdSQS. Both 5-epi-4′-O-demethylancistrobertsonine C and ancistroealaine B formed one hydrogen bond each through the Val139 residue. Ancistrobenomine A interacted with Val139 and Ser13 residues through H-bonds. The maximum number of H-bonds (5 H-bonds) was formed by E5700 with the residues Cys255, Gln259, Gln176 and Val139 within the binding pocket of LdSQS. Three intermolecular H-bonds were observed for 3-O-methyldiplacol and the residues involved were found to be Tyr36, Gln259 and Cys255.

H-bond, carbon hydrogen bond, hydrophobic and electrostatic interactions were the most significant interactions identified between top hit ligands and the LdSQS protein. The hydrophobic interactions were largely based on Pi-Sigma, Pi-Pi stacked, Pi-Pi T shaped, Alkyl and Pi-Alkyl types. Mainly, Cys255, Gln259, Glu176, Val139, Ser13 and Tyr36 were attributed in H-bond formation in protein-ligand complexes. All ligands were associated with hydrophobic interactions especially via Val143, Phe16, Tyr36, Ala140, Leu175, Leu147, Pro258, Leu39, Val139, Phe35, Phe254, Met171, Met114, Val19, Met118, Cys255 whereas Val43, Leu175 and Ala140 residues involved consistently in forming hydrophobic bonds with each top hit ligand.

Details of specific intermolecular interactions between eight ligands studied and LdSQS protein are depicted in Table S3. The 3D interaction diagrams and 2D interaction plots of the best binding conformation from each ligand (top hits) are illustrated in Figure 5. The representation of top hit ligands occupy the same binding cavity is shown in Figure 5I.

Figure 5

The best docking pose of top hit compounds in the binding cavity of LdSQS. A) Ancistrotanzanine B, B) E5700, C) Ancistrotanzanine A, D) Ancistrobenomine A,E) 5-epi-4′-O-demethylancistrobertsonine C, F) Ancistroealaine B, G) Ramiflorine B, H) 3-O-methyldiplacol. The 3D interaction diagrams are shown at the left whereas 2D interaction plots are depicted at the right. I) The superimposition of top hit ligands on each other. The ligands, Ancistrotanzanine A (purple), Ramiflorine B (yellow), E5700 (cyan), Ancistrobenomine A (green), 5-epi-4′-O-demethylancistrobertsonine C (orange), Ancistroealaine B (gray), 3-O-methyldiplacol (magenta) and Ancistrotanzanine B (white) occupy in the same binding pocket of LdSQS protein. Green sphere represents the Magnesium ion.

Pharmacophore modeling was performed using PharmaGist server in order to elucidate the 3D pharmacophoric features of each top hit ligand. These chemical features are necessary to interact with a specific biological target structure either to trigger or to block its biological response. Top hit ligands used in this study were bound at the same binding pocket in the same orientation. Therefore analysis of pharmacophores of each ligand is useful to identify lead compounds to discover better inhibitors against LdSQS [70]. The recognized features of pharmacophores (chemical features) for each top hit ligand are summarized in Table S4. The highest number of pharmacophores was found to be 27 for 3-O-methyldiplacol while the least number was recorded as 15 for reported inhibitor, E5700.

The pairwise structural alignment details (Table 3) demonstrated common chemical features shared between E5700 (pivot molecule) and each top hit phytochemical. Ancistrotanzanine A shared the maximum number of common features (7) with E5700 although the best docked phytochemical; ancistrotanzanine B shared the least number of features (3). The pairwise structural alignment of E5700 with each phytochemical (top hits) are depicted in Figure 6.

Table 3

Pairwise structural alignment showing common pharmacophoric features between E5700 (pivot molecule) and top hit compounds.

Score	F	SF	Ar	Hydrophobic	Donors	Acceptors	N	P	Molecules
9.33099	6	5	1	1	1	3	0	0	E5700-Ancistrobenomine A
9.32368	5	4	2	1	1	1	0	0	E5700-Ancistroealaine B
9.01744	4	3	2	0	1	1	0	0	E5700-5-epi-4′-O-demethyl ancistrobertsonine C
8.43524	7	7	1	3	1	2	0	0	E5700-Ancistrotanzanine A
8.134	5	5	2	2	0	1	0	0	E5700-Ramiflorine B
7.52356	4	3	1	0	1	2	0	0	E5700-3-O-methyldiplacol
7.51348	3	3	2	0	0	1	0	0	E5700-Ancistrotanzanine B

F- Features, SF- Spatial features, Ar- Aromatic, N- Negatives, P- Positives.

Figure 6

Structural alignment of pivot molecule, E5700 (gray) with (A) Ancistrotanzanine A (red), (B) Ancistroealaine B (magenta), (C) 3-O-methyldiplacol (violet), (D) Ancistrotanzanine B (blue), (E) Ancistrobenomine A (cyan), (F) 5-epi-4′-O-demethylancistrobertsonine C (green) and (G) Ramiflorine B (yellow). E5700 is displayed in ball and stick style whereas phytochemicals are shown in sticks style. White spheres represent hydrogen bond donors, yellow spheres represent hydrogen bond acceptors, and magenta spheres specify aromatic features and green spheres depict hydrophobic features. B) & F) structures, hydrogen bond acceptor and hydrogen bond donor are in the same position. C) & E) structures, hydrogen bond donor and hydrogen bond acceptor are in the same position.

Pharmacokinetic parameters of drugs are strongly influenced by their physicochemical properties. Consequently, selected compounds were evaluated for their physicochemical properties through Molinspiration tool. They are listed in Table 4.

Table 4

Physicochemical properties and drug-likeness of top hit compounds.

Description	miLogP	TPSA (A°)	nNHA	MW (g/mol)	nHBA	nHBD	nRotB
Ancistrotanzanine B	5.86	49.30	31	419.52	5	0	5
E5700∗	2.99	69.06	32	433.55	6	2	4
Ancistrotanzanine A	5.56	60.29	30	405.49	5	1	4
Ancistrobenomine A	5.08	81.05	31	419.48	6	2	5
5-epi-4′-O-demethyl ancistrobertsonine C	5.67	51.17	31	421.54	5	1	4
Ancistroealaine B	5.42	59.96	30	407.51	5	2	4
Ramiflorine B	5.79	56.08	35	466.63	5	3	3
3-O-methyldiplacol	4.91	116.45	33	454.52	7	4	7

miLogP- octanol/water partition-coefficient, TPSA- Topological Polar Surface Area, nNHA-number of Non Hydrogen Atoms, MW- Molecular Weight, nHBA-number of Hydrogen Bond Acceptors, nHBD-number of Hydrogen Bond Donors, nRotB- number of Rotatable Bonds, Volume- Volume of molecule, Lipinski's rule- (MW =<500, HBA =<10, miLogP =<5, HBD =<5), Veber rule- (nRotB =<10, TPSA =<140).

E5700∗- Reported inhibitor 1.

One of the criteria to determine drug oral bioavailability is the molecule to satisfy the Lipinski's rule of 5 (LRo5) with no more than one violation. The reported inhibitor and 3-O-methyldiplacol satisfied all LRo5 and Veber rules. The following phytochemicals viz. ancistrotanzanine A, ancistrotanzanine B, ancistrobenomine A, 5-epi-4′-O-demethylancistrobertsonine C, ancistroealaine B and ramiflorine B showed high logP value than 5, violating only one criterion of LRo5. But all compounds theoretically may demonstrate good oral bioavailability because none of them violated more than one criterion in LRo5. There are several orally administered medicines that are in clinical practice, for treatment of cancer, hepatitis C and cardiovascular diseases which violate the LRo5 [71].

Early investigation of in silico pharmacokinetic and toxicity properties is essential to avoid the chance of drug failures during clinical trials. Hence each top hit ligand was evaluated to predict the pharmacokinetic descriptors and toxicity profile using the pkCSM tool (Table 5).

Table 5

In silico ADMET predictions of top hit compounds computed by pkCSM server.

Description	Ancistrotanzanine B	E5700	Ancistrotanzanine A	Ancistrobenomine A
Absorption
HIA∗ (%)	98.598	97.491	96.695	96.275
Water solubility [log (mol/L)]	-6.354	-2.66	-5.733	-5.256
Distribution
VDss∗∗ [log (L/kg)]	0.441	0.959	0.219	-0.562
BBB#permeability (logBB)	-0.527	-0.51	-0.446	-0.666
Metabolism
CYP2D61 substrate	No	No	No	No
CYP3A41 substrate	Yes	Yes	Yes	Yes
CYP1A21 inhibitor	Yes	No	Yes	Yes
CYP2C191 inhibitor	Yes	No	Yes	Yes
CYP2C91 inhibitor	Yes	No	Yes	Yes
CYP2D61 inhibitor	No	No	No	No
CYP3A41 inhibitor	Yes	Yes	Yes	Yes
Excretion
Total clearance [log (ml/min/kg)]	0.505	0.702	0.521	0.406
Renal OCT2## substrate	No	Yes	No	No
Toxicity
AMES toxicity	No	No	No	Yes
hERG∗∗∗I inhibitor	No	No	No	No
hERG∗∗∗ II inhibitor	Yes	Yes	Yes	Yes
Hepatotoxicity	No	Yes	No	Yes

HIA∗- Human Intestinal Absorption, VDss∗∗- steady state Volume of Distribution, BBB#- Blood-Brain Barrier, CYP1- Cytochrome P450 isoenzymes, OCT2##- Organic Cation Transporter 2, hERG∗∗∗- human ether-a-go-go gene.

Accepted range of values.

HIA %- >30%.

Water solubility [log (mol/L)] - insoluble < -10 < poorly soluble < -6

VDss [log (L/kg)] - high>0.45, low < -0.15.

BBB permeability (logBB) - poorly permeable to the brain < -1, readily cross the blood-brain barrier>0.3.

Total Clearance [log (ml/min/kg)]- <1.28.

All selected ligands possessed high human intestinal absorption percentages ranging from 87.185 % to 98.598% demonstrating that they can be easily absorbed by the intestine. Water solubility predictions are given in log (mol/L) (insoluble < -10 < poorly soluble < -672]. Total clearance and substrate of renal organic cation transporter 2(OCT2), are the most important descriptors to detect excretion of drug molecules. All selected molecules exhibited acceptable values for total clearance with respect to standard value (<1.28). Total clearance is measured as a combination of renal clearance and hepatic clearance. There is an inverse relationship between total drug clearance and drug half-life where the higher the total clearance value, the lower is the half-life of drug. Therefore the drug molecules with higher clearance can be eliminated with no significant action on biological targets. The organic cation transporter 2 transports positively charged drugs from the blood supplying the proximal tubule to the inside of the tubular epithelial cell. As a result, molecules that are substrate of OCT2 can be eliminated faster than non-substrates. Two compounds namely 5-epi-4′-O-demethylancistrobertsonine C and reported inhibitor were found to be substrates for OCT2 among other top hits. Finally toxicity analysis was done to check Ames toxicity, cardiotoxicity and hepatotoxicity potential of each selected compound. Out of eight top hits, three were positive for Ames toxicity that may lead to DNA mutations whereas five compounds showed no potential for hepatotoxicity. Interestingly all compounds were negative for hERG1 channel inhibition indicating cardio protective nature of the top hits.

Experimental results from previous studies have reported E5700 and ER119884 as potent SQS inhibitors [9, 10]. Two compounds had shown antileishmanial activity against promastigotes and intracellular amastigotes from different incubation periods of L. amazonensis Josefa. These two compounds were used as references in virtual screening. The IC50 values of E5700 were reported as 14.7 nM (promastigotes) and 4 nM (amastigotes) whereas ER119884 revealed IC50 values of 1.7 nM (promastigotes) and 0.9 nM (amastigotes) for 72 h incubated samples [10]. After a 96 h incubation, E5700 showed IC50 values of 6.7 nM (promastigotes) and 2.0 nM (amastigotes) whereas ER119884 displayed IC50 values of 0.9nM (promastigotes) and 0.5 nM (amastigotes) [10].

According to the published data, ancistrotanzanine B, ancistrotanzanine A, ancistrobenomine A, 5-epi-4′-O-demethylancistrobertsonine C and ancistroealaine B showed antileishmanial activity against L. donovani parasite in terms of IC50 values of 1.6 μg/ml (=3.81 μM), 1.8 μg/ml (=4.44 μM), >10 μg/ml (≥23.84 μM), 12.1 μg/ml (=28.70 μM) and 10 μg/ml (=24.54 μM) [32, 33, 34, 35]. Antitrypanosomal activity of aforementioned compounds found to be IC50 values of 1.5 μg/ml (=3.58 μM), 1.7 μg/ml (=4.19 μM), 4.8 μg/ml (=11.44 μM), 4.7 μg/ml (=11.15 μM) and 17.6 μg/ml (=43.19 μM) [32, 33, 34, 35]. In addition to that ramiflorine B exhibited median lethal dose (LD50) value of 12.63 ± 5.52 μg/ml (=27.07 μM) against L. amazonensis promastigotes [41]. Moreover, 3-O-methyldiplacol inhibited proliferation of axenic L. donovani and T. brucei brucei cell lines at an IC50 of 7.2 ± 1.1 μg/ml (=15.84 ± 2.42 μM) and 7.2 ± 1.7 μg/ml (=15.84 ± 3.74 μM) respectively [42].

Previous reports showed different susceptibilities of miltefosine in terms of IC50 values at the micromolar range against intracellular amastigotes of five Leishmania spp. viz. L. amazonensis (3.21 μM), L. braziliensis (5.40 μM), L. guyanensis (4.02 μM), L. chagasi (4.46 μM) and L. donovani (0.22 μM) [73]. Another study has shown IC50 values of miltefosine ranging from 3.74-6.15 μM against L. donovani promastigotes isolated from Indian VL patients [74].

Conclusion

The docking results revealed seven phytochemicals surpassed the threshold value of -9.00 kcal/mol and strongly bound to the LdSQS. The highest docking score was observed for ancistrotanzanine B (-9.83 kcal/mol) whereas 3-O-methyldiplacol formed three hydrogen bonds, binding tightly to the active site of LdSQS.

Through multiple sequence alignment, two aspartate rich regions and two conserved motifs were identified among protein sequences of Leishmania species. These motifs contain the binding site and catalytic site residues of the protein (Figure 1B). The LdSQS sequence showed close evolutionary relationship with L. mexicana.

PharmaGist server was applied to find pharmacophores of each top hit compound and common chemical features shared between each phytochemical with E5700. The phytochemical, 3-O-methyldiplacol fell within the permissible range of all the ADMET and physicochemical parameters compared to the reference compounds. These top hit phytochemicals can be further structurally optimized and test in vitro and in vivo concentrations to be developed as antileishmanial drugs.

Declarations

Author contribution statement

Padmika Madushanka Wadanambi: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Uthpali Mannapperuma: Analyzed and interpreted the data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the commercial, public, or not-for-profit sectors.

Data availability statement

Data included in article/supp. material/referenced in article.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.