Bulbs of E. plicata were collected in Tracuateua, Pará, Brazil (Lat. 1.1436°, Long. 46.9551°), a specimen was deposited at Museu Paraense Emílio Goeldi (MG 202631). The research project complies with national guidelines and international legislation, registered on the platform of the National System of Management and Genetic Heritage and Associated Traditional Knowledge (SISGEN), under registration number A49DEEE that grants license for collecting the species. The EEEp was obtained by macerating the dry powder of the bulbs (924 g) in ethanol (2 L for 7 days), subjected to fractionation in an open chromatographic column, using silica gel mesh (63–200 mm) as stationary phase and mobile solvents of increasing polarity (hexane, dichloromethane, ethyl acetate and methanol), obtaining the fractions: FHEp, FDEp, FAEEp, and FMEp, concentrated in a rotary evaporator. The FDEp was subjected to fractionation in thin layer chromatography on a preparative scale and showed four yellow spots with different retention factors, which were removed separately, being named Subfraction FA1, FA2, FA3, and FA4. From FA3, isoeleutherin was isolated, and identified by nuclear magnetic resonance (NMR) spectra, using a Bruker Advance DPX 400 MHz NMR spectrometer (Bruker Ascend).

Inhibition of Leishmania growth was evaluated in vitro by cultivating promastigotes of L. amazonensis in stationary phase (5×10⁶ parasites) in the presence of the extract, fractions and isoeleutherin (200–3.125 μg/mL) in 96-well culture plates (Nunc, Nunclon^®, Roskilde, Denmark), for 72 h at 26°C. Viability was assessed by measuring the cleavage of MTT [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl tetrazolium bromide] (Sigma) (Ngure et al., 2009). Absorbances were measured using a multiwell scanning spectrophotometer (Molecular Devices, Spectra Max Plus, Canada) at 490 nm. Amphotericin B (AmpB) was used as a positive control (25–0.3906 μg/mL). The concentration of products required to inhibit 50% of the viability of L. amazonensis (IC₅₀) was determined by applying a sigmoidal regression of the individual concentration response curves of the compounds. The data are representative of independent experiments, carried out in triplicate, which presented similar results (Mota et al., 2011).

Murine macrophages (4 × 10⁵ cells) were seeded on round glass coverslips into 24-well culture plates (Nunc) in RPMI 1640 medium (Sigma and Aldrich), supplemented with 20% fetal bovine serum (FBS), 2 mM L-glutamine, penicillin 50 IU/mL, and streptomycin 50 μg/mL, pH 7.4. After 24 h of incubation at 35°C in 5% CO₂, promastigotes of L. amazonensis in stationary phase were added to the wells (4 × 10⁶ parasites) to promote infection of macrophages, and the cultures were incubated for 4 h at 35°C in 5% CO₂. Next, the free parasites were removed by extensive washing with RPMI 1640 medium, and the infected macrophages were quantified and treated with the extract, FA2, FA3, fraction FA4, and isoeleutherin (500–125 μg/mL, each) for 72 h at 35°C in 5% CO₂. The negative control consisted of infected macrophages and culture medium. The positive control used AmpB (100–25 μg/mL). Then, coverslips were stained with Giemsa, and the percentage of inhibition of intramacrophage viability of Leishmania was determined by counting the number of amastigotes per 100 macrophages on each coverslip, under a light microscope (×100 magnification). The presented data were performed in triplicate, and the IC₅₀ was determined using GraphPad Prism version 5.04 (Silva, 2005).

Free online platform PreADMET (https://preadmet.webservice.bmdrc.org/) was used to predict absorption, distribution, metabolism, and elimination properties. Molecular descriptors related to Lipinski’s rule of five (Lipinski et al., 2001) and Veber extensions (Veber et al., 2002), such as molar mass (MW), octanol/water partition coefficient (cLog P), number of hydrogen bond acceptor (HBA) and hydrogen bond donor (HBD), number of rotatable bonds (nRotb) and topological polar surface area (TPSA) were predicted in Molinspiration (https://www.molinspiration.com) (Molinspiration, 2023).

Bioavailability prediction considered Lipinski’s “Rule of Five”, where a good drug candidate will have molecular weight <500, partition coefficient (log P) < 5, no more than five hydrogen bond donors and ten acceptors of hydrogen bonding (Lipinski et al., 2001). In the pharmacokinetic studies, intestinal absorption (Human Intestinal Absorption = HIA) was evaluated, considering parameters in the range of 0%–20% (low absorption), 20%–70% (moderate absorption), >70% (high absorption; (Hou et al., 2007). Molecule permeability in Caco-2 and MDCK cells were considered high when presented values >70 nm/s, average of 4–70 nm/s and low <4 nm/s (Yee, 1997; Yazdanian et al., 1998; Balimane et al., 2000). For distribution analysis, values >90% indicate strong binding to albumin, while <90% binding will be moderate to weak (Sun et al., 2018). As for the ability to cross the blood-brain barrier, the following criteria were used: freely crosses BBB >2.0, moderately crosses BBB values between 2.0 and 0.1 and reduced crossing or not crosses <0.1 (Ajay and Murcko, 1999).

Molecular docking was used to explore the possible conformations of the ligand with the binding receptor, estimating the intensity of enzyme-ligand interaction (Meng et al., 2011). Isoeleutherin and the subset derived from isoeleutherin analogue compounds were prepared from SMILES files downloaded from the ZINC (Irwin et al., 2012) and eMolecules (www.emolecules.com) databases (eMolecules, 2023). The crystallographic structure of the TR enzyme was retrieved from the Protein Data Bank (PDB) under the code 2JK6 (Baiocco et al., 2009) with a resolution of 2.95Å, prepared using the Chimera program (Lang et al., 2009), removing water molecules, ligands and adding hydrogen atoms. Molecular docking simulations were performed using the GOLD 2020.1 program (Cambridge Crystallo-graphic Data Center—CCDC, Cambridge, United Kingdom), which uses a genetic algorithm to generate and select conformations of flexible compounds that bind to the receptor site of a protein (Jones et al., 1995).

Compounds were scored by applying the GoldScore scoring function with a 100% efficient search. The binding site was defined based on studies by Padey et al. (2016), at a 10Å sphere centered on the flavin-adenine dinucleotide ligand (FAD); (Pandey et al., 2017a; Pandey et al., 2017b). The methodology was validated through redocking, evaluated using the fconv 1.24 program (Neudert and Klebe, 2011). The procedure was carried out to evaluate the convergence of the results and to determine the smallest value of the root mean square deviation (RMSD), making it possible to select the best interactions of the crystallographic ligand in the complex formed. To analyze the hydrogen bonds and hydrophobic interactions between the selected ligands and the enzyme’s amino acids, the PoseView online server was used (Stierand et al., 2006), a tool that displays molecular complexes that incorporate a simple and easy-to-perceive arrangement of the interactions formed between ligands and amino acids (Stierand and Rarey, 2007). Quinacrine was used as a control drug, which was chosen because it is a small molecule like the ones studied and because it has activity on TR (Saravanamuthu et al., 2004).

Firstly, the calculation of electrostatic potential charge was performed for the structure obtained during docking using the Gaussian 03 program (Frisch et al., 2003), applying the Restricted Electrostatic Potential (RESP) (Bayly et al., 1993) associated with the Hartree-Fock methods (Fock, 1930) and HF/6-31G (d,p) base function (Hariharan and Pople, 1973). The protonation states of all amino acid residues present in the enzyme were determined at pH 7.0 using the Propka server (Olsson et al., 2011). Molecular dynamics simulations were conducted using the AMBER 18 program, implemented by the pmemd. CUDA module (Li et al., 2013). The general AMBER force field (GAFF) (Wang et al., 2004) was employed to describe ligands, and the MMFF99SB force field (Hornak et al., 2006) was used for enzyme amino acid residues. The enzyme-ligand complex was solvated with the TIP3P explicit solvent model in a cubic box with an edge length of 12 Å, with the inclusion of Clˉ counter-ions to achieve electrical neutrality in the system, using the tleap module included in the AMBER program (Jorgensen et al., 1983).

Minimizations and heating were carried out with a SANDER module (Case et al., 2005), the systems were fragmented into four stages of energy minimization. In the first phase, 25,000 minimization steps were performed, divided into 10,000 steps, performed with the steepest descent method and 15,000 with the conjugate gradient. The remaining three phases were performed in 10,000 minimization steps, using the same methodology for each step. Subsequently, the systems were gradually heated using the Langevin algorithm, with protein atoms subjected to a restriction constant of 25 kcal/mol. Å2, considering the NVT set from 0 K to 298 K (Loncharich et al., 1992).

Periodic boundary conditions were simulated using the Particle Mesh Ewald (PME) method, employed for long-range electrostatic interactions (Darden et al., 1993). Cutoff distances for the long-range and van der Waals interactions were set at 9 Å. After minimization and equilibrium of the system, during the period of 50 ns of DM simulation, an integration time of 2.0 fs was produced using the Verlet algorithm (Verlet, 1968), considering the NPT set adjusted to a temperature of 298 K and pressure of 1 atm for each enzymatic binding complexes, all bonds with hydrogen atoms were restricted using the SHAKE algorithm (Ryckaert et al., 1977). The CPPTRAJ module of AMBERTOOLS 18 (Roe and Cheatham, 2013) was used to carry out the structural analyzes of RMSD and B-factor. The RMSD calculation verified the stability of the systems in relation to the initial structure, and the application of the B-factor identified which amino acid residues were most flexible in the enzyme. The binding free energy calculation was performed in the last 10 ns, using the AMBERTOOLS 18 modules CPPTRAJ and MMPBSA. py (Miller et al., 2012). The single-path MM-PB(GB)/SA protocol considers identical conformation states of the protein-ligand complex, unbound protein and free ligand (Kollman et al., 2000; Massova and Kollman, 2000).

The binding free energy (ΔG_bind) of the complexes was determined according to Eq. 1, where ΔH is the enthalpy term, TΔS represents the product between absolute temperature and entropy resulting from the conformations obtained in the DM simulation, ΔEMM consists of the energy obtained by molecular mechanics and ΔGbind, solv is the free energy of solvation. Entropy contributions to energy result from changes in translation, rotation, and vibration. ∆ G b i n d   = ∆ H − T ∆ S ≈ ∆ E M M   + ∆ G b i n d , s o l v   − T ∆ S (1)

The mechanical energies, represented by Eq. 2, are calculated involving the individual contributions of internal energy (ΔEint), electrostatic (ΔEele) and van der Waals (ΔEvdW). ∆ E M M   = ∆ E int   + ∆ E e l e   + ∆ E v d w (2)

The sum to compose the internal energy, described in Eq. 3, contains the contributions of bond length (ΔEbond), bond angles (ΔEangle) and torsion angles (ΔEtorsion). ∆ E int   = ∆ E b o n d   + ∆ E a n g l e   + ∆ E t o r s i o n (3)

Solvation free energy (ΔG_{bonding, solv}) results from the sum of polar (ΔG_PB/GB) and non-polar (ΔG_n-polar) contributions. According to Eq. 4, the polar electrostatic contribution to the solvation free energy can be calculated by the Poisson-Boltzmann method (PB) or by the generalized methods of Born approximation (GB). ΔG_PB and ΔG_GB were calculated using the generalized Born model (igb = 2) (Onufriev et al., 2000) and in the MMPBSA method, considering the dielectric constants of solute (Ajay and Murcko, 1999) and solvent (Yazdanian et al., 1998; Sun et al., 2014). ∆ G b i n d , s o l v   = ∆ G P B / G B   + ∆ G n o n − p o l a r (4)

The nonpolar energy is estimated by the product of the surface tension γ with a value equal to 0.0072 kcal/mol Ų and the surface accessible solvent area (SASA) according to Eq. 5. ∆ G n − p o l a r = γ ∙ S A S A (5)

The interaction free energy that allows determining the individual contribution of all residues to the free energy of the complexes, obtained by the MMGBSA method, is described in Eq. 6, elucidating its importance in the active site of the enzyme (Alves et al., 2020). ∆ G l i g a n d − r e s i d u e = ∆ G v d W + ∆ G e l e + ∆ G G B + ∆ G S A (6)

The main residues that contribute to the total energy were visualized using the CHEWD plugin (Raza et al., 2019) of the Chimera program (Pettersen et al., 2004). The van der Waals (ΔG_vdW) and electrostatic (ΔG_ele) interactions between amino acid residues of the TR enzyme were determined by the SANDER module, implemented in the AMBER 18 program.

From the ethanolic extract (EEEp; yield 2.84%) the following fractions were obtained: hexane (FHEp; yield 3.1%), dichloromethane (FDEp; yield 19.6%), ethyl acetate (FAEEp; yield 10.2%) and methanol (FMEp; 61.9% yield). Isoeleutherin was isolated from FDEp and its identification is described by Borges et al. (2020).

In vitro anti-leishmanial assays were conducted to evaluate the effect of the samples (EEEp, FDEp, FAEEp, FMEp, and isoeleutherin) on promastigote and intracellular amastigote forms of L. amazonensis. The results indicated no activity for the extract and its fractions, with an inhibitory concentration of 50% (IC₅₀) exceeding 200 μg/mL, while isoeleutherin demonstrated activity (IC₅₀ = 25 μg/mL). Additionally, it is important to note that the concentrations of the samples used in this study were not cytotoxic to macrophages (CC₅₀ > 500 μg/mL), similar to the control drug (amphotericin B; CC₅₀ > 100 μg/mL). The EEEp, its fractions, and isoeleutherin showed CC₅₀ greater than 500 μg/mL and IC₅₀ greater than 200 μg/mL against L. amazonensis, which is why the selectivity index calculation was not possible.

Then, macrophages infected with promastigote form were submitted to treatment with EEEp, FDEp, FAEEp, FMEp, and isoeleutherin (concentration of 500, 250 e 125 μg/mL). In Figure 2, amastigotes forms can be observed around the destroyed cells, exposed to EEEp (Figure 2A) and its fractions FDEp (Figure 2B), FAEEp (Figure 2C), FMEp (Figure 2D) and isoeleutherin (Figure 2E). The destruction of macrophages observed, likely caused by intracellular forms of L. amazonensis, since EEEp, its fractions, and isoeleutherin did not show cytotoxicity on macrophages. No reductions in the number of amastigotes were observed in infected cells when compared to the negative control (Figure 3).

Amastigote forms can be observed around the destroyed cells exposed to treatment with EEEp, FDEp, FAEp, FMEp, and isoeleutherin, highlighting that fractionation did not enhance the activity against amastigote forms of L. amazonensis, as the fractions and isoeleutherin remained inactive. The large number of promastigote forms around the destroyed macrophages indicates that they were unable to invade the cells, considering that these forms can survive in unfavorable conditions. Some of them assumed a more rounded shape, resembling amastigotes (Figure 2).

Isoeleutherin (ISO; Figure 4) was employed as the starting molecule in the search for analogs through structural similarity in databases, with the aim of examining how structural modifications would impact physicochemical and pharmacokinetic properties. As a result of this process, 20 analogs (Figure 4) were identified and subsequently subjected to in silico prediction studies.

Isoeleutherin and its 20 analogues did not violate Lipinski’s rule, showing octanol-water partition coefficient (miLog p ≤ 5); Molecular Mass (MM) ≤ 500 g/mol, number of hydrogen bond acceptor groups (nHAG) ≤ 10 and number of hydrogen bond donor groups (nHDG ≤5) (Lipinski et al., 1997). The values obtained for the topological polar surface area (TPSA) parameters ≤140Ų; number of rotational bonds (Nrotb) ≤ 10 and molecular volume are within the limits established by Veber’s descriptors (Table 1) (Veber et al., 2002).

The results suggest isoeleutherin and its analogues have moderate permeability in Caco2 and low to moderate in MDCK, whereas the permeability of isoeleutherin in MDCK is moderate. However, it is possible to observe that these compounds have high intestinal absorption (AIH). Isoeleutherine appears to bind to albumin moderately, with some structural alterations increasing the affinity for albumin. Only isoeleutherin and compound 18 freely crossed the BBB, structural alterations reduced the ability to cross the barrier. All compounds appear to be metabolized by CYP3A4, with variations in the magnitude of metabolism. Structural changes did not significantly impact the CYP inhibitory potential (Table 2).

Therefore, it is important to evaluate the impact of structural alterations in different aspects, such as: physicochemical, pharmacokinetic and receptor binding. Among the isoeleutherin analogues, there were no significant alterations in physicochemical and pharmacokinetic aspects, except in distribution, where alterations in plasma protein binding and distribution to the CNS were observed.

Molecular docking of the isoeleutherin and its 20 analogues were evaluated at the catalytic site of the TR enzyme (PDB 2JK6), performed at a distance of 10Å from the Flavin-Adenine Dinucleotide (FAD) cofactor, with RMSD values below 2Å, demonstrating that the redocking was successful according to literature data (Hevener et al., 2009). Connections that occurred in the active site of the enzyme were analyzed using the PoseView online server (Stierand et al., 2006), taking into account the interactions performed with residues Cys52, Cys57, His461’, and Glu466’, involved in the redox metabolism of Leishmania (Baiocco et al., 2009). The GoldScore scoring function was employed to predict the binding affinity of the most stable configuration, and the highest scores were considered to select the top 10 compounds with potential anti-leishmanial activity, as shown in Table 3.

Quinacrine is an effective and widely used antiparasitic drug with potential adverse effect. In this experiment, it was employed as the control drug. Compound 17 (Zinc317780204) was selected among the isoeleutherin analogues, based on the results found in molecular docking, showing that the molecule’s conformation interacts with the enzyme, through hydrophobic interactions with the amino acid residues Thr51 and hydrogen bonds with the Tyr198, Cys57, Lys60 residue. Isoeleuterin presented conformations with the enzyme by hydrophobic interactions with the residue Gly161 and hydrogen bonds with residues Ser14, Thr51, Ser162, and Arg287 (Figure 5).

The coordinates obtained from the molecular docking simulation procedures, representing the most stable configurations of quinacrine, isoeleutherin, and the analogue compound (Zinc317780204), were subjected to 50 ns of Molecular Dynamics (MD) to assess the structural stability of the enzyme-ligand complexes through the analysis of RMSD values over time, as depicted in Figure 6. After the process was initiated, the RMSD values of the TR-ligand complexes fluctuated (1–4 Å) until reaching equilibrium (10 ns) and then remained stable throughout the simulation, with average values of 2.50 and a standard deviation of ±0.31 for quinacrine, 2.70 Å and a standard deviation of ±0.39 for isoeleutherin, and 2.33 Å and a standard deviation of ±0.23 for Zinc317780204. This indicates that the compounds remained within the enzyme cavity, showing minimal conformational changes in the complex structure.

The flexibility of the protein regions in relation to each ligand was evaluated through the graph of B-factor (Figure 7). The greatest fluctuations occurred in regions corresponding to the bands of amino acid residues Asn91-Gly80 (highlighted in red), Glu410-Thr397 (highlighted in yellow) and Ser489 fragment region, highlighted in green, considered the most flexible regions of the enzyme. These residues are found close to the Cys52 and Cys57 active sites and play an important functional role in enzyme inhibition (Verma et al., 2012; Pandey et al., 2016).

Comparing the binding free energy values of the two methods, described in Table 4, both are favorable and stable for complex formation. The calculation of binding free energy allowed for quantifying the affinity between the enzyme-ligand system. When comparing the free energy values obtained through the MM-GB(PB)SA method, it was observed that quinacrine, used as the reference drug in the present study, displayed the lowest value. The compound Zinc317780204 followed, with values lower than those of isoeleutherin, indicating that these compounds may be useful in situations where the reference drug has limitations.

Energy decomposition by residue was employed to identify the TR enzyme residues involved in interactions with Zinc317780204, isoeleutherin, and quinacrine, based on the energy contribution calculated by the MMGBSA method for the last 10 ns of MD simulations. Free binding energy values below −1 kcal/mol were set as the criteria for selecting amino acid residues that participated in the most favorable interactions, contributing to ligand stabilization within the complex.

The analysis of residue decomposition values reveals those that played a significant role in the total interaction energy and stabilization of the TR-Zinc317780204 (Cys52, Tyr198, Val332, and Thr335), TR-isoeleutherin (Cys52, 57, and Thr335), and TR-Quinacrine (Cys52, Tyr51, Cys57, Arg287, Asp327, Met333, and Ala338) complexes, contributing significantly to the overall free energy of the complex (Figure 8).