Machine Learning and Structure-based Virtual Screening for Selecting Cinnamic Acid Derivatives as Leishmania major DHFR-TS Inhibitors

1. Introduction

Leishmaniasis is a neglected tropical disease (NTD) caused by protozoan parasites of the genus Leishmania, which are transmitted by the bite of infected sandflies. This disease affects millions of people worldwide, particularly in developing countries with poor health infrastructure. The primary clinical forms of the disease are visceral, cutaneous, and mucocutaneous. According to the World Health Organization (WHO), the global burden of leishmaniasis is estimated to be around 700,000 to 1 million new cases each year, with 90% of the cases occurring in just six countries: Afghanistan, Algeria, Brazil, Colombia, Iran, and Syria [1,2]. The sandflies that transmit leishmaniasis are most active at night and breed in wet soil, organic matter, or animal burrows [3]. In Colombia, 10 out of the 20 species that can infect both humans and other living beings are present. The cutaneous leishmaniasis (CL) form is the most frequent (98-99%), with the population under five years old and immunocompromised individuals being the most affected [4,5]. The number of CL cases reported in Colombia in 2022 was 4,906, with the departments of Amazonas, Boyacá, Caquetá, Cesar, Córdoba, Cundinamarca, Putumayo, Santander, and Sucre being the most affected areas [6].

Since the late 1980s, Leishmania-HIV co-infection has been reported in 35 countries, and there have also been other cases of Leishmania-Malaria co-infection, which are associated with the worsening of the clinical condition of patients with leishmaniasis. This co-infection type has increased the disease's burden due to the greater difficulty of clinical treatment [7,8]. Currently, antimonial compounds are the primary treatment for leishmaniasis; however, they present high toxicity and resistance in some endemic regions. To address these challenges, alternative drugs have been developed, such as liposomal amphotericin B, which significantly reduces the side effects and treatment duration associated with free amphotericin B, but is expensive [9,10]. Other drugs, such as paromomycin and miltefosine, have been associated with high toxicity, resistance, and teratogenic and abortive effects, promoting the discovery and development of low-cost, highly effective drugs with low toxicity [11]. Furthermore, it is worth noting that while Leishmania is a parasitic disease mainly affecting humans, it also affects animals such as dogs and rodents, which can serve as reservoirs for the parasite and increase the risk of transmission to humans [12,13]. Therefore, efforts to develop effective treatments and control measures must be considered.

High-throughput screening (HTS) has been used since the early 1990s to test the activity of large numbers of molecules against different diseases and thereby identify potential hits for drug development [14]. However, the uncertainty of success, as well as the time and screening costs, limit the use of this technique [15]. In recent years, chemoinformatics tools (e.g., molecular docking, machine learning) have been utilized to conduct in-silico studies that can predict the interactions between a protein and a ligand, reducing the number of actual laboratory experiments and accelerating the drug discovery process more efficiently and cost-effectively [14, 16]. The different research conducted in this field has led to the development of increasingly efficient and better-classifying models, which take advantage of large compound databases, opening the possibility of studying diseases that mainly affect poorer populations (NTD), which are not attractive to large industries and big pharma [17].

Leishmaniasis is commonly treated with plants from the Asteraceae family in traditional medicine. Given the diversity of this family (32,913 species) and the wide range of phytochemicals they contain, including alkaloids, coumarins, flavonoids, benzofurans, sterols, and terpenoids, they are considered a promising source of new leishmanicidal compounds [18]. Some secondary metabolites studied in this family have been sesquiterpenoids [19,20], triterpenes [21], phytosterols [22], and kauranes [23]. However, although they have shown activity to inhibit the disease, their pIC₅₀ is not large enough, and compounds that are effective at low concentrations and selective against the parasite are preferred. A group of compounds that have not yet been studied, with records reporting promising in-vitro activity, are derivatives of cinnamic acid belonging to the Asteraceae family [24-26].

Gouri et al. report some natural inhibitors against Leishmania amastigotes such as luteolin (IC₅₀ = 3.12 μM), quercetin (IC₅₀ = 10.5 μM), chrysin (IC₅₀ = 13 μM), apigenin, myricetin, cinnamic acid (IC₅₀ = 0.25 μM), licochalcone A (IC₅₀ = 0.9 μM), which can play an important role in drug discovery [24]. Peixoto et al., on the other hand, evaluated the biological activity of 25 cinnamic acid derivatives against Leishmania braziliensis amastigotes, obtaining promising results and finding that aromatic rings with an oxygen as a heteroatom had a beneficial effect in terms of activity against Leishmania [25]. Considering that heterocyclic compounds have been of great importance for drug development in the pharmaceutical industry, derivatives of cinnamic acid, which is an aromatic carboxylic acid commonly substituted in trans position by an acrylic acid group, represent an interesting starting point for directing studies in the search for possible hits against different species of leishmaniasis [27]. Although some of these compounds have already been studied, many more remain to be analyzed.

In the present study, a computational approach was carried out to select potential inhibitors of the bifunctional enzyme dihydrofolate reductase-thymidylate synthase (DHFR-TS) of Leishmania major given that it plays a crucial role in the synthesis of DNA in trypanosomatids, which is essential for the parasite's reproduction [28]. In this regard, a custom-made, in-house library containing 314 specialized metabolites derived from cinnamic acid was virtually screened. Initially, a ligand-based predictive classification model was developed using experimental information on the IC₅₀ values retrieved from in-vitro assays of reported compounds against Leishmania. In parallel, using a hybrid LmDHFR-TS model constructed based on its amino acid sequence [29], a structure-based ranking, through molecular docking calculations, was performed using the investigated specialized metabolite database. Through a consensus analysis, molecules with the highest probability of being inhibitors by both approaches were classified as possible hits. These secondary metabolites were evaluated through in-vitro assays using the recombinant LmDHFR-TS, and ADMET properties were calculated to determine their pharmacokinetic properties.

2. Results and Discussion

2.1. Combined Ligand-/Structure-Based Virtual Screening Approach using LmDHFR-TS

2.1.1. Ligand-Based Virtual Screening

Initially, a bank of compounds with inhibitory activity reported against LmDHFR-TS was created from the ChEMBL database. These molecules were classified as active or inactive based on their reported IC₅₀ values, using a cut-off point of IC₅₀ = 5.0. Duplicate molecules were removed to perform the data curation process and obtain a virtual screening model with high prediction efficiency. Additionally, molecules with an IC₅₀ value within ±0.1 of the cut-off point were considered. Finally, a total of 790 molecules were selected for model training. Among these, 378 were inactive (47.8%) and 412 were active (52.2%).

In the ligand-based process, VolSurf+ (128) and AlvaDesc (more than 4000) molecular descriptors were calculated from the three-dimensional representation of each database compound. These descriptors were then used to create the Random Forest (RF) model in Knime software (KNIME 4.5.0, the Konstanz Information Miner, Copyright 2003–2014, www.knime.org), consisting of 200 decision trees. The Gini Index was employed as a split criterion in the RF model to reduce the number of false-positive results. The dataset underwent a five-fold cross-validation strategy, where it was divided into five subsets, each containing an 80% modeling set and a 20% validation set. The modeling set was exclusively used for model construction and further subdivided into multiple training and test sets, maintaining an 80%/20% split ratio. These procedures were conducted following the approach described by Fourches et al. [30].

Molecular descriptors play a crucial role in drug discovery and development, as they represent the molecular and chemical properties of the studied compounds. In this study, the employed descriptors were instrumental. VolSurf+ generates three-dimensional (3D) molecular descriptors based on the distribution of molecular electrostatic potentials and hydrophobicity, encompassing molecular surface properties such as size, shape, and electrostatic potential distribution [31,32]. On the other hand, AlvaDesc offers a wide range of descriptor types, including constitutional descriptors (the number and type of atoms, bonds, and functional groups in the molecule), topological descriptors (representing molecular shape, size, and complexity), electrostatic descriptors (conveying molecular polarity and charge distribution), and quantum mechanical descriptors (relating to the electronic structure and properties of the molecule) [33,34].

The performance of the RF model was evaluated to compare the effectiveness of the two types of descriptors. This evaluation involved calculating classification precision, recall, F1-score, and Matthew's correlation coefficient (MCC). Additionally, Receiver Operating Characteristic (ROC) curves were analyzed, and the Area Under the ROC Curve (AUC) was calculated (see Figure 1). These evaluation metrics are commonly used to measure the performance of binary classification models. ROC curves and their AUC are often employed to assess the performance of models that produce continuous output scores or probabilities. AUC serves as a scalar measure of the model's overall ability to distinguish between positive and negative cases [33,35]

According to the parameters displayed in Figure 1, it is evident that the MCC and AUC values for both the test sets and cross-validation are higher for AlvaDesc descriptors compared to those obtained for VolSurf descriptors. However, considering that a higher AUC value indicates a more remarkable classification ability of the model and that MCC is expressed in a range of -1 to 1 (where a high value close to 1 suggests a strong correlation between the predicted class and the true class), good values were obtained for both AlvaDesc (AUC: 0.863 and 0.906, MCC: 0.554 and 0.645) and VolSurf (AUC: 0.855 and 0.884, MCC: 0.539 and 0.598) descriptors.

Regarding precision, recall, and F1-score, good and similar values were obtained for both models, except for the recall for inactive compounds in the model created using VolSurf descriptors, which was low with a value of 0.69. Sensitivity and specificity measures were also calculated to assess the performance of the RF model. For AlvaDesc, the values were 0.807 and 0.752, while for VolSurf, the values were 0.843 and 0.690, respectively. These results indicate a tendency to have few false negatives, a higher value of true negatives, and a lower false-positive rate for both descriptors. The applicability domain was also determined, confirming that the regression model can provide reliable predictions.

Ligand-based virtual screening (VS) was then employed to predict the potential inhibitory activity of 314 compounds derived from cinnamic acid in the Asteraceae family, as documented in the literature. Figure 2 displays the structure and probability of the five best compounds classified using AlvaDesc descriptors. These compounds were (E)-2-hydroxy-3',6'-dimethoxychalcone (103) [36], apigenin 7-O-(6´´-caffeoyl)-glucoside (235) [37], montamine (63) [38], 3-O-p-coumaroyl-betulinic acid (150) [39], and cordoin (202) [40]. In addition, Figure 2 presents the top five compounds predicted using VolSurf descriptors: 6,8-di-C-β-glucopyranosylchrysin (242) [41], montamine (63) [38], dihydrocubebin (305) [42], prebalanophonin (312) [43], and 4-O-feruloyl 5-O-caffeoylquinic acid (96) [44].

Among all the tested compounds, 116 were classified as active using AlvaDesc molecular descriptors, with probability values ranging from 0.50 to 0.71. On the other hand, 93 compounds were considered active with VolSurf molecular descriptors, and their probability values ranged from 0.50 to 0.86. Some of these molecules were previously reported to exhibit various activities, such as analgesic activity (305), antimalarial activity (150), cytotoxic activity (63), acting as anticancer agents (202), and demonstrating antiproliferative properties (312) [38,40-42,45-47].

Regarding the best compounds, only one contains nitrogen in its structure (63). The rest have various oxygen atoms, forming heterocycles or containing carbonyl groups, ethers, and alcohols. Additionally, one of them is a steroid (150), and another is glycosylated (242).

2.1.2. Structure-Based Virtual Screening

The Structure-Based Virtual Screening (VS) was conducted utilizing a hybrid homology model of LmDHFR-TS [29], a bifunctional enzyme with a critical role in the metabolic pathway of Leishmania parasites, as well as several protozoa species. Leihsmania genus is autotrophic for folate and unconjugated pteridines, with the enzyme DHFR-TS playing a pivotal role in the reduction of dihydrofolate to tetrahydrofolate, a cofactor in the biosynthesis of thymine in nucleotide metabolism [48].

To evaluate the putative inhibitory capability of cinnamic acid derivatives against LmDHFR-TS, molecular docking calculations were performed using Molegro software. The results were validated by redocking the co-crystalized ligand, i.e., ethyl 4-(5-{[(2,4-diaminoquinazolin-6-yl)methyl]amino}-2-methoxyphenoxy)butanoate (DQ1) along with the reference inhibitor methotrexate (MTX) (Figure 3).

The compounds were ranked according to the predicted docking binding energy using the probability calculation shown below (Eq. 1), as previously reported by Herrera-Acevedo et al. [19, 20]. The ten compounds exhibiting the highest probability of being active are presented in Table 1. Ranked compounds that did not previously show high ligand-based probability values but appeared among the best-ranked derivatives through a structure-based approximation are represented in Figure 4, along with their respective structure-based probability (P_SB) values.

$$P_{SB}=(E_i/E_{min})>0.5\mathrm{and}{E}_i<E_{ligand}$$

(1)

Where $P_{SB}$ is the structure-based probability; $E_i$ is the docking energy of compound $i$, where $i$ ranges from 1 to 314 (cinnamic acid derivatives dataset); $E_{min}$ is the lowest energy value of the dataset; and $E_{ligand}$ is the ligand energy from co-crystalized inhibitor.

The results showed that energy-based scoring values were lower for the cinnamic acid derivatives than for the reference ligands, suggesting that the studied compounds have a better affinity with the LmDHFR-TS active site in the molecular recognition process. Additionally, the docking results showed that 24.5% of the 314 cinnamic acid derivatives dataset had $P_{SB}$ values above 0.5, and 64 of these top-ranked compounds had a lower docking score than methotrexate, which obtained -114.15 kJ/mol.

Three of the top-ranked molecules predicted to have high ligand-based probability values based on the RF model also demonstrated high structure-based probability values. Specifically, compound 242, ranked fourth in the structure-based classification (Table 1), was the best classified in the ligand-based VS model with VolSurf descriptors. Compounds 235 and 63, positioned among the top ten compounds in structure-based VS with docking scores of -161.4 kJ/mol and -160.1 kJ/mol, respectively, also showed high ligand-based probabilities. Compound 235 was predicted to be the second-best structure with high potential for inhibition using the model built with AlvaDesc descriptors, while Compound 63 was classified in the top three for both RF models (AlvaDesc and VolSurf molecular descriptors).

The analysis of residues for the best poses in the top three compounds showed that the residues responsible for ligand binding (Val30, Val31, Ala32, Ile45, Trp47, Asp52, Met53, Phe56, Val87, Pro88, Fhe91, Leu94, Val156, Tyr162, and Thr180) have been previously reported in the literature as part of the active site [49]. Certain characteristics of these residues, such as accessibility and charge distribution, enable selective drug design against these protozoans without affecting human enzymes [49]. The interaction diagrams in Figure 5 illustrate that the compound with the highest docking score (compound 241, Figure 5C) possesses heterocyclic rings similar to the reference ligands, with oxygen atoms replacing the nitrogen atoms present in the reference ligands. However, due to the similar electronegativities of nitrogen and oxygen, these atoms favor nearly identical interactions with the enzyme's active site.

Compounds 164 and 21 lack heterocyclic rings but contain benzene rings, which participate in π–π and π–alkyl interactions. Additionally, these compounds exhibit a relevant number of oxygen-containing groups, such as esters, ethers, and carboxylic acids, facilitating interactions with both residues within the active site and other residues. Specifically, the carboxylic moiety facilitates Van Der Waals interactions, crucial as they occur with the amino groups in the reference ligands and appear to be important since they are present in the three top-ranked molecules. On the other hand, compounds 242 and 140, containing only hydroxyl groups, are less favorable in this binding mode. Although both compounds are isomeric, compound 164 has few favorable interactions (8 interactions), and compound 21 has more interactions (25 interactions), but it also has some unfavorable interactions (Ala32, Ala47, Ala119).

Furthermore, it was observed that all ligands adopted a U-shaped conformation like the reference ligands DQ1 and MTX (Figure 5F), and most of them have strong hydrogen bonding interactions with the enzyme (Val156, Val30, Lys95, Met53, Phe91, and Arg97), which are important determinants for binding [50]. To further explore this behavior, a topological polar surface area (TPSA) map was constructed for both the reference ligands and the best-ranked compounds (Figure 6).

The results of the TPSA maps confirmed a similar spatial distribution among the three top-ranked compounds with respect DQ1 and MTX. An electron-deficient region was identified at the top of the molecule (Figure 6, blue area), which is consistently present in all evaluated molecules, including the two reference ligands. This observation rationalized the similar binding behavior within the active site of LmDHFR-TS, particularly with Met53 as a common crucial contact for these test compounds.

2.1.3. Consensus Analysis of the Two VS Approaches

A combined approach was employed to ascertain the potential activity of cinnamic acids against the LmDHFR-TS enzyme and to mitigate the selection of false-positive compounds. This approach incorporated probability scores derived from both structure-based and ligand-based virtual screening (VS) methods in conjunction with the true-negative rate obtained from the RF model (Eq. 2) [19].

The design of this approach aimed to assign a higher weight to the ligand-based probability scores (considering their reliance on experimental pIC₅₀ values), in contrast to the structure-based probability scores, which are founded on protein-ligand interactions. This weighting scheme significantly reduces the risk of incorrectly classifying inactive molecules as active (false positives) [23].

$${CA}_{Lm}=\frac{\left[P_{SB}+\left(\left(1+{TN}_{LB\left(AD\right)}\right)x{P}_{LB\left(AD\right)}+\left(1+{TN}_{LB\left(VS\right)}\right)x{P}_{LB\left(VS\right)}\right)\right]}{\left[3+{TN}_{LB\left(AD\right)}+{TN}_{LB\left(VS\right)}\right]}$$

(2)

Where ${CA}_{Lm}$= combined-approach probability, $P_{SB}$= structure-based probability, TN = true-negative rate, and $P_{LB}$= ligand-based probability (AD=AlvaDesc descriptors and VS=VolSurf descriptors).

Table 2 presents the results of the best-ranked compounds calculated from the consensus analysis equation. The compounds ranked among the top five for each method are highlighted in bold. Except for 235, all compounds were classified as potentially active in all virtual screening approximations used in this study. The consensus analysis identified 110 compounds with combined-approach probability values greater than 0.5; however, only 47% of these compounds (52) were classified as active through the three in-silico models used in this study (Table S1). Compound 63 (montamine) was the top-ranked compound. Montamine is an indole alkaloid that has been isolated from Asteraceae species, such as Centaurea schischkinii and Centaurea montana. Previous studies have reported its anticancer properties [38,51], but its efficacy against Leishmania has not been investigated.

The second best-ranked compound is 6,8-di-C-β-glucopyranosylchrysin (242), a derivative of chrysin obtained from Lychnophora ericoides (Asteraceae). Compared to compound 69 (chrysin) and 231 (techtochrysin), classified as inactive, the glycosylated derivative 242 has more hydroxyl groups, enabling interactions with the enzyme's active site. In previous studies, chrysin was biofunctionalized with gold particles due to its low bioavailability, poor absorption, and rapid excretion issues, aiming to neutralize Leishmania parasites through its activity against the kinase-3 enzyme [52]. However, compound 242 could represent an alternative due to its hydrophilic character resulting from the glycosyl groups, potentially inhibiting Leishmania parasites by interacting with LmDHFR-TS.

The third- and fourth best-ranked compounds were 4-O-feruloyl-5-O-caffeoylquinic acid (96) and lucenin-2, 6,8-di-C-β-glucopyranosylluteolin (241), both extracted from the genus Lychnophora, specifically Lychnophora ericoides [41] and Lychnophora salicifolia [44], respectively. Additionally, apigenin 7-O-rutinoside (39), lithospermic acid (237), diarctigenin (306), and isolappaol A (308) – four cinnamic acid derivatives that previously exhibited moderate values in both RF models and the molecular docking calculations (all classified as active) – appeared among the top ten ranked compounds in the combined approach (Figure 7). Hence, these compounds emerge as interesting antileishmanial candidates, as they exhibit activity across all models and maintain consistency in their probability values. Notably, consensus scoring methods are known to enhance hit rates by diminishing the likelihood of false positives [53].

The compounds 4-(3,4-dihydroxybenzyl)-2-(3,4-dihydroxyphenyl)tetrahydrofuran-3-carboxy-O-β-D-glucopyranoside (306) and 7-(3,4-dihydroxyphenyl)-3′,4′-dihydroxy-7,8,7′,8′-tetrahydronaphtho[8,8′-c]furan-1(3H)-one (308) are two lignans found in certain species of Asteraceae. Notably, Hypochaeris radicata (native to Europe, northern Asia, and parts of North Africa) and Arctium lappa (native to Europe and Asia) have been reported as a natural sources of these compounds. However, A. lappa is widely disseminated in America and H. radicata has also become invasive in regions as far-flung as New Zealand and Chile [54]. Conversely, compound 237, lithospermic acid, is a common polycyclic phenolic carboxylic acid that has been isolated from species of multiple botanical families, including Lamiaceae and Asteraceae. It has demonstrated a wide range of beneficial properties, acting against cardiovascular diseases and hepatitis. It allows endothelium-dependent vasodilatation, lowers blood pressure, and produces antioxidant effects [55,56].

2.2. Molecular Dynamics Simulations

Molecular dynamics (MD) studies were conducted to assess the protein–ligand stabilities, considering factors such as solvent, ions, pressure, and temperature for compounds 237, 306, and 308. These three compounds were identified as potential inhibitors of LmDHFR-TS through the consensus analysis of the approaches outlined in this study. Methotrexate (MTX) was also used as the reference ligand.

The structural stability was evaluated through root-mean-square deviation (RMSD) measurements. During the simulated time (100 ns), all tested compounds showed similar behavior concerning the apoenzyme of LmDHFR-TS (apoLmDHFR-TS, the protein devoid of a ligand) and the LmDHFR-TS···MTX complex. After 20 ns, all ligands (including the reference ligand) displayed diminished disturbances, reaching values ranging from 0.23 to 0.27 nm. The derivative 237 exhibited a higher disturbance between 60 and 85 ns, with RMSD values varying close to 0.1 nm (Figure 8a).

Scheme 8. b), allowing us to examine residue flexibility in the presence of various ligands [23]. The RMSF plot revealed higher fluctuations in the N-terminus region, reaching approximately 0.35 nm, while the C-terminal residues exhibited lower fluctuations with values close to 0.20 nm. Regions characterized by well-defined tertiary structures, such as α-helices or β-sheets, consistently displayed RMSF values ranging from 0.10 to 0.20 nm.

In general terms, all evaluated compounds exhibited the same behavior; however, we identified some specific cases. The residues Glu218 and Thr410, located in the protein's loop regions, showed the highest fluctuations for the apoenzyme, with Glu218 having approximately twice the RMSF value compared to the complexes with MTX and the tested cinnamic acid derivatives. Among the selected compounds, compound 237 displayed higher fluctuations in the loop regions than the other derivatives and MTX, with Gly118, Arg254, and Arg380 being the most variable amino acids. Compounds 306 and 308 exhibited a similar behavior throughout the simulation, with reduced flexibility when complexed with LmDHFR-TS.

The structural compactness and level of mobility of the protein-ligand complexes throughout the simulation period were evaluated using the radius of gyration (RoG) plot (Figure 8c) [23]. In the case of LmDHFR-TS, during the initial half of the 50 ns simulation, the complexes with cinnamic acid derivatives displayed RoG values that were indistinguishable from those of the control MTX and apoLmDHFR-TS (ranging from 2.64 nm to 2.70 nm). This finding suggests a high level of stability and low fluctuations in the tertiary structure. However, after 60 ns, compounds 237, 306, and 308 exhibited a similar behavior (varying between 2.64 nm to 2.70 nm) with increased perturbation compared to the DHFR-TS···MTX complex and the apoenzyme, which maintained a consistent mean value with fluctuations ranging from 2.62 to 2.64 nm.

After conducting molecular dynamic simulations, binding free energies for complexes involving compounds 237, 306, 308, and the control (MTX) with LmDHFR-TS were determined using the MM/PBSA method. The complexes of both benzylbutyrolactone-type lignans (306 and 308) and the polyphenolic acid (compound 237) with LmDHFR-TS showed binding free energies of -111.1 kJ/mol, -81.0 kJ/mol, and -91.6 kJ/mol, respectively. In all cases, the energy was higher than the -124.5 kJ/mol observed for the complex of MTX with LmDHFR-TS (Table 3).

All evaluated complexes, including the reference MTX, showed the same contribution pattern with negative energy values from van der Waals, electrostatic, and solvent-accessible surface area (SASA) parameters to the binding free energy. The van der Waals parameter presented the highest negative contribution, with values lower than -209 kJ/mol. This observation suggested that non-polar electrostatic interactions are the primary factor responsible for the molecular recognition of the LmDHFR-TS binding site by test compounds.

Regarding polar solvation, all compounds positively contributed to the total binding energy, with similar values for compounds 237, 308, and MTX. Conversely, diarctigenin (306) contributed less to this parameter. Moreover, electrostatic interactions negatively contributed to the binding free energies, being more relevant for MTX, with a negative contribution of -57.3 kJ/mol. Meanwhile, its influence was between 35% and 50% for the evaluated cinnamic acid derivatives compared to the reference ligand.

2.3. In-vitro Enzymatic Activity Inhibition for selected cinnamic acid derivatives (compounds 237, 306, and 308) against LmDHFR-TS and HsDHFR.

Five compounds from our in-house compound library were available to validate the results derived from the combined approach employing the two virtual screening (VS) methodologies through in-vitro enzymatic inhibition assays. Compounds 237, 306, and 308, classified as active in all approaches, along with hesperidin (140), an interesting flavonoid previously reported to exhibit antileishmanial activity through the induction of apoptosis and inhibition of sterol C-24 reductase [57], were included in this in-vitro exploration. Isovitexin 4′-O-glucoside and rutin, which displayed a moderate level of activity and were classified as inactive in one of the three approaches, were also evaluated against LmDHFR-TS, with methotrexate as the positive control.

The IC₅₀ values were determined by analyzing concentration-response curves in the 0.1–128 μM range, using spectrophotometric monitoring of enzymatic activity under a standard DHFR assay. This analysis yielded a spectrum of values spanning from 6.1 to 53.2 μM, corresponding to pIC₅₀ values ranging from 4.27 to 5.21. The compounds 237, 306, and 308 were the most active molecules against LmDHFR-TS. Hesperidin (IC₅₀ = 21.6 μM) showed high activity against the target among the three tested flavonoids, with observed IC₅₀ values of 53.2 μM and 41.7 μM for isovitexin 4′-O-glucoside and rutin, respectively (Table 4).

Structurally, we correlated the inhibitory activity against LmDHFR-TS with the relation between H-bond acceptors and donors (carbonyl and hydroxyl groups). Among the lignans 306 and 308, the presence of the γ-butyrolactone moiety showed that the most active compound (306) has a higher number of carbonyl groups compared to 308—a feature also observed in Lithospermic acid (237). Nevertheless, the low inhibitory activities of the glycosylated flavonoids (hesperidin, isovitexin 4′-O-glucoside, and rutin) indicate that the high presence of hydroxyl groups seems to have a negative effect on the inhibitory activity.

Subsequently, based on the obtained results from in-vitro tests using the recombinant protein of Homo sapiens (Hs) DHFR, the selectivity index (SI) was calculated. The IC₅₀ values obtained against HsDHFR showed a different pattern, suggesting different mechanisms of action for these two proteins. Moderate SI values were revealed, with both benzylbutyrolactone-type lignans (compounds 306 and 308) showing the highest SI values at 4.6 and 4.4, respectively. Both lignans exhibited higher SI values than MTX, which was used as a positive control (Table 4).

2.4. Pharmacokinetic properties predictions.

Pharmacokinetic properties, i.e., absorption, distribution, metabolism, excretion, and toxicity (ADMET) of the compounds 237, 306, and 308, were predicted by using ADMETlab 2.0 and OSIRIS Datawarrior 5.5.0 [58,59]. Different approaches were used to assess the likelihood of oral bioavailability, and the results were mixed: on the one hand, all compounds were compliant with Lipinski's "rule of five" [60], but on the other hand, none of the above mentioned were accepted under the parameters established by Pfizer [61] and GSK [62]. Consequently, these compounds might not perform well regarding oral bioavailability (Table S2).

Concerning cytochrome P450 (CYP) and its isoenzymes, compound 237 was considerably likely to inhibit CYP2C9, and so were compounds 306 and 308 with CYP2C19, CYP2C9, and CYP3A4; this means that the metabolism of other drugs might be affected by the compounds evaluated herein. Meanwhile, compound 237 is likely to behave as a substrate for CYP2C9, and compounds 306 and 308 with CYP1A2, CYP2C19, CYP2C9, CYP2D6, CYP3A4; thus, these compounds can be metabolized by the mentioned isoenzymes. Additionally, none of the studied compounds possess mutagenic, tumorigenic, reproductive, or irritant effects. Identifying potential hERG channel blockers has great importance when assessing the risk of cardiotoxicity [63]; for the three structures, the probabilities of hERG blocking were 0.212 at the most.

3. Materials and Methods

3.1. Cinnamic acid derivatives in-house dataset

A custom-made, in-house virtual library of 314 distinct cinnamic acid derivatives was built from 76 scientific articles using various search criteria, including keywords such as Asteraceae, Cinnamic Acid Derivatives, Lignans, Polyphenols, Flavonoids, and others. ChemAxon MarvinSketch [ChemAxon, version 21.18.0 (2021), a calculation module developed by ChemAxon, https://www.chemaxon.com/, accessed on 12 January 2023] was used to design all the structures.

The three-dimensional (3D) structures for the entire set were generated using Standardizer software [JChem, version 21.18.0 (2021), a calculation module developed by ChemAxon, https://www.chemaxon.com/, accessed on 12 January 2023]. This software standardized the structures, added hydrogens, performed aromatic form conversions, and refined molecular graphs in three dimensions. The final dataset was saved in Special Data File (.sdf) format.

3.2. Classificatory Machine Learning Models.

The analyses described below utilized Knime 4.5.0 software (KNIME 4.5.0, the Konstanz Information Miner, Copyright 2003–2014, www.knime.org) [64]. The process commenced with importing those descriptors generated by the Volsurf+ [31,32] and AlvaDesc [33,34] programs in CSV format.

Subsequently, these descriptors underwent segmentation via the 'Partitioning' node, implementing the stratified sampling option, with 80% of the initial dataset designated as the training set and the remaining 20% composing the test set. Random splits were also explored while maintaining consistent ratios for both training and test sets.

The model's creation processes entailed utilizing the modeling set and the RF algorithm, executed through a five-fold cross-validation procedure employing WEKA nodes. This approach provides a robust and efficient means to evaluate a model's performance by partitioning the data into five subsets for testing and training, facilitating model selection and generalization assessment [23].

RF models were configured with 200 trees, a random number of generator seed of 1, and employed the Gini Index as the split criterion for both training and cross-validation sets. Performance analysis of the selected models involved an assessment of their internal and external aspects, incorporating parameters such as sensitivity (true-positive rate), specificity (true-negative rate), and accuracy (overall predictability), derived from the confusion matrix.

To provide a more comprehensive understanding of the model's performance beyond accuracy, the ROC curve was employed, generated through a 'ROC curve' node, which relies on sensitivity and specificity. AUC values were obtained from this curve, where a value of 0.5 indicates an inability to distinguish between the two groups, and an AUC of 1 signifies perfect separation without overlap [65]. Furthermore, the MCC was calculated, with a value of 1 representing perfect prediction, 0 denoting random prediction, and -1 indicating complete disagreement between prediction and observation [66].

Additionally, a performance evaluation of the RF model using AlvaDesc and VolfSurf+ descriptors was conducted, incorporating precision, recall, and F1-score for both active and inactive sets.

3.3. Molecular docking calculations.

Molecular docking calculations involved the hybrid model of LmDHFR-TS bound to methotrexate (PDB ID: MTX) [29] and the three-dimensional structures of the cinnamic acid derivatives. We conducted these calculations using Molegro 6.0.1 software.

To ensure consistency, we removed all water molecules from both the enzyme and compound structures, and we prepared them to use the software's default settings. The MolDock scoring function was utilized, considering internal ES, internal H-bond, and Sp2-Sp2 torsions as criteria for evaluating the ligands.

The molecular docking process was executed through 10 runs utilizing the MolDock SE algorithm. It allowed for a maximum of 1500 interactions, maintained a population size of 50, included up to 300 steps, employed a neighbor distance factor of 1.00, and returned a maximum of 5 poses. To cover the enzyme's ligand-binding site, we established a grid with a 15 Å radius and 0.30 Å resolution [23,29].

Our results were categorized according to docking scores, reflecting the free energy or affinity of the interactions. Each calculation was repeated three times to ensure reliability. For comparison, we employed methotrexate (MTX) as a control.

Total Polar Surface Area (TPSA) maps were calculated using Spartan 14 for Windows Spartan'14 (Wavefunction Inc., Irvine, CA, USA) [67], and visualization of two-dimensional residual interaction diagrams was accomplished using Discovery Studio Visualizer v21.1.0.20298 (BIOVIA, Dassault Systèmes, San Diego, CA, USA) [23,29].

3.4. Molecular dynamics simulations.

Molecular dynamics (MD) simulations were conducted in YASARA Structure v. 19.12.14 [68], employing the AMBER14 force field to model the enzyme and ligand–enzyme systems. Before the simulations, each protein underwent hydrogen bond optimization, and chloride (Cl⁻) and (Na⁺) ions were added to the model systems through the transferable intermolecular potential 3-point (TIP3P) employing 0.997 g/L density for solvating the simulation cell. Acid dissociation constant values (pKa) were calculated for enzyme’s titratable amino acids with H-bonding network and the side-chain placement with a rotamer library (SCWRL) algorithm. Periodic boundary conditions were applied to facilitate the simulations, involving a cell size set 10 Å larger than the protein's size in all instances. An initial 5000-cycle energy minimization step was carried out using the steepest gradient approach. MD simulations used the particle-mesh Ewald (PME) method to account for long-range electrostatic interactions (8-Å cut-off distance). The simulations were performed under physiological conditions at 298 °K, pH 7.4, and 0.9% NaCl. Temperature control was maintained using a Berendsen thermostat while keeping the pressure constant. A multiple-time step algorithm with a time step of 2.00 fs was employed. Finally, MD simulations were run for 100 ns under constant pressure and the Berendsen thermostat, with snapshots saved at intervals of 100 p, using the YASARA macro (md_run.mcr) for all simulation phases. Subsequent analyses were also carried out using the default YASARA macro scripts. The molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method was employed to calculate the binding free energies of apoenzyme and enzyme–ligand complexes from the resulting MD trajectories using the g_mmpbsa tool in Gromacs 5.0.5 (open source, http://www.gromacs.org) [69] on an Ubuntu 12.04 server, using NPT and periodic boundary conditions, as previously reported [29, 70].

3.5. LmDHFR-TS and HsDHFR enzymatic inhibition assays.

Purification and kinetic characterization of the recombinant LmDHFR-TS protein were performed according to the previously reported procedures [29,71-72], while HsDHFR protein was obtained from the commercial assay kit (CS0340, Merck KGaA, Darmstadt, Germany). Thus, the in-vitro evaluation of top-ranked selected compounds (237, 306, 308, hesperidin, rutin, and isovitexin 4′-O-glucoside) for inhibitory activity against LmDHFR-TS and HsDHFR was conducted using a spectrophotometric assay under standard DHFR conditions. These tested compounds were available from our in-house compound library. Rutin, lithospermic acid, and rutin were commercially purchased (>98%, Merck KGaA, Darmstadt, Germany). Isolappaol and diarctigenin were isolated from a commercial A. lappa powdered root extract (Prescribed For Life, Fredericksburg, TX, USA) through successive column chromatography, whose spectroscopic data was identical to those of previous reports [73-74].

The assay consisted of LmDHFR-TS or HsDHFR (2.7 nM), bovine serum albumin (BSA, 1 mg/mL), N-[tris(hydroxymethyl)-methyl]-2-aminoethanesulfonic acid (TES) buffer (100 mM, pH 7.0, 150 mM β-mercaptoethanol, 2 mM ethylenediaminetetraacetic acid (EDTA)), and nicotinamide adenine dinucleotide phosphate (NADPH, 100 μM) with varying concentrations of the test compounds (0.1-128 μM). The reaction was initiated by adding the substrate (7,8-dihydrofolate (H2F), 20 μM) and was monitored for 360 seconds at 340 nm (i.e., oxidation of NADPH to NADP⁺) to determine the initial reaction rate (Vo) through linear regression analysis of the resulting absorbance profile. All measurements were performed in triplicate, and MTX was used as a positive control [29]. The resulting Vo values were used to calculate the % inhibition, as 100 - (Ri / Rc x 100), where Ri is the Vo in the presence of the inhibitor, and Rc is the Vo in the absence of inhibitors (1% DMSO v/v final concentration). The % inhibition was measured for at least five concentrations (0.1-128 μM) for each test compound (i.e., cinnamic acid derivatives and MTX), and concentration-response curves (% inhibition vs. Log[inhibitor]) were constructed using non-linear regression to determine the IC₅₀ using GraphPad Prism 7.0 (GraphPad, San Diego, CA, USA) [29].

3.6. Pharmacokinetic properties predictions.

For the compounds 237, 306, and 308, the ADMET parameters were calculated using the ADMETlab 2.0, an integrated online platform for predictions of ADMET properties [58]. Drug toxicity prediction was performed using OSIRIS DataWarrior v.5.2.1, based on the following parameters: mutagenicity, tumorigenicity, reproductive effect, and irritability [59].

4. Conclusions

This study identified three cinnamic acid derivatives, Lithospermic acid (237), Diarctigenin (306), and Isolappaol A (308), as potential inhibitors of LmDHFR-TS using a combined Virtual Screening approach (structure/ligand-based). Two random forest models were built using different molecular descriptors. Sensitivity and specificity measures were obtained to evaluate the RF model performance. The models classified 116 (AlvaDesc) and 93 compounds (VolSurf) as active, showing a tendency to minimize false negatives.

Molecular docking revealed that 24.5% of the 314 cinnamic acid derivatives had values above 0.5, with 64 of them having a lower docking score than methotrexate, the reference ligand. A consensus analysis combining the RF models with molecular docking identified 110 compounds with combined-approach probability values greater than 0.5. From them, 47% were classified as active through the in-silico models, identifying some compounds with potential leishmanicidal activity that a single approach had not previously highlighted. Lithospermic acid (237), diarctigenin (306), and isolappaol A (308) were among the top-ranked compounds, and their binding mode was evaluated using molecular dynamics. Finally, in-vitro assays using recombinant LmDHFR-TS validated the computational results, with 237, 306, and 308 exhibiting significant activity against LmDHFR-TS. However, moderate selective indices (SIs) were observed when assays were performed using HsDHFR. Despite this finding, higher SI values than MTX were observed. Thus, these three tested compounds emerged as an interesting alternative as hits against LmDHFR-TS; however, specific assays against the parasitic forms of Leishmania major are required to extend a clearer prospect for fighting this neglected tropical disease.