Alkyl and Aryl Derivatives Based on p-Coumaric Acid Modification and Inhibitory Action against Leishmania braziliensis and Plasmodium falciparum

In low-income populations, neglected diseases are the principal cause of mortality. Of these, leishmaniasis and malaria, being parasitic, protozoan infections, affect millions of people worldwide and are creating a public health problem. The present work evaluates the leishmanicidal and antiplasmodial action of a series of twelve p-coumaric acid derivatives. Of the tested derivatives, eight presented antiparasitic activities 1–3, 8–12. The hexyl p-coumarate derivative (9) (4.14 ± 0.55 μg/mL; selectivity index (SI) = 2.72) showed the highest leishmanicidal potency against the Leishmania braziliensis amastigote form. The results of the molecular docking study suggest that this compound inhibits aldehyde dehydrogenase (ALDH), mitogen-activated kinase protein (MPK4), and DNA topoisomerase 2 (TOP2), all of which are key enzymes in the development of Leishmania braziliensis. The data indicate that these enzymes interact via Van der Waals bonds, hydrophobic interactions, and hydrogen bonds with phenolic and aliphatic parts of this same compound. Of the other compounds analyzed, methyl p-coumarate (64.59 ± 2.89 μg/mL; IS = 0.1) demonstrated bioactivity against Plasmodium falciparum. The study reveals that esters presenting a p-coumarate substructure are promising for use in synthesis of derivatives with good antiparasitic profiles.


Introduction
Neglected diseases typically affect tropical and subtropical regions, principally developing countries [1,2], with poor and marginalized people [3] being the most affected. According to the World Health Organization, neglected diseases affect more than 1 billion people around the world [4]. These diseases cause great impacts on public health and include Chagas disease, leishmaniasis, malaria, schistosomiasis, onchocerciasis, lymphatic filariasis, trypanosomiasis and dengue [5].

Antileishmanial and Antiplasmodial Activity of Compounds 1-12
The results were evaluated using the 50% effective concentration (EC 50 ) for parasitic activity. We analyzed the collection (twelve compounds) for activity against Leishmania braziliensis in its amastigote form, against Plasmodium falciparum, and for cytotoxicity against U-937 macrophages using 50% lethal concentration (LC 50 ). Of the derivatives analyzed in Table 1, eight compounds (1-3, and 8-12) showed moderate to good activity against Leishmania braziliensis amastigotes. However, the twelve compounds exhibited weak antiplasmodial activity. Table 1. Antileishmanial and antiplasmodial activity, cytotoxicity and molar mass of compounds 1-12.
Compounds presenting three side chain carbons (3, 4, and 5) exhibited differing EC 50 values (respectively 6.3 ± 1.14 µg/mL, 28.33 ± 0.27 µg/mL, and 86.98 ± 18.13 µg/mL). This may be due to an alteration in the n-propyl radical (3) to isopropyl (4), contributing to the decrease in leishmanicidal activity seen in further branching in compound 4 [47,48]. Compound 3 was fourteen times more active than compound 5. The presence of the oxygen atom in the compound 5 side chain likely contributes to a reduction in antiparasitic activity. Meanwhile, Nóbrega et al. (2019) [45] observed that the presence of side chain oxygen in a methoxy-ethyl derivative (3,4,5-trimethoxycinnamate) increases leishmanicidal activity. The leishmanicidal activity values of both compounds 4 and 5 were lower (more potent) than their cytotoxic activity values against U-937 macrophage cells (respectively 96.23 ± 18.89 µg/mL, 126.71 ± 16.5 µg/mL. These same cytotoxicity values were also much higher (less potent) than that of compound 3 (11.89 ± 0.8 µg/mL) which had high cytotoxicity.
Hexyl p-coumarate (9), dodecyl p-coumarate (10), 4-methylbenzyl p-coumarate (11) and 4-isopropylbenzyl p-coumarate (12), (respectively 4.14 ± 0.55 µg/mL, 16.34 ± 0.49 µg/mL, 5.69 ± 0.3 µg/mL and 11.21 ± 0.2 µg/mL) displayed high leishmanicidal activity, likely due to their relatively high lipophilicity. Yet Otero et al. (2014) [52] found in their study that hexyl cinnamate presents less leishmanicidal activity than dodecyl cinnamate. However, in the present study 9 was considered the most potent compound compared to both compound 10 and the other derivatives analyzed. The results obtained are in accordance with the preparation plan for the compounds, in which alkyl and benzyl substituents were added to increase the lipophilicity of the compounds. In this way, compounds with good absorption were obtained to cross the parasite membrane resulting in greater bioavailability [53][54][55]. Several mechanisms of action can occur, such as greater access to organelles that are biological targets for leishmanicidal activity [56]. As suggested in the molecular docking study, hexyl p-coumarate should inhibit the enzymes ALDH, MPK4, and TOP2 (see Table 2) which are essential for Leishmania braziliensis growth. Regarding the substitution of the aliphatic chain in the para position of the aromatic ring in compounds 11 and 12, a difference in leishmanicidal activity was noted since compound 11 was twice as potent as compound 12 which possessed the more bulky substituent. These same results were also observed by Lopes et al. (2019) [44] in relation to trypanocidal activity.
If the side chain of the p-coumarate derivative is increased, an increase in the lipophilicity of the compounds is observed for pentyl p-coumarate of (7), isopentyl p-coumarate of (8), hexyl p-coumarate of (9) and dodecyl p-coumarate (10). As for antiplasmodial activity, there was a difference in EC 50 between these compounds; compound 8 (95.48 ± 1.84 µg/mL) was more potent than compound 7 (156.37 ± 8.2 µg/mL). Yet Arsianti et al. (2017) [59] observed that alkyl ester groups presenting linear chains such as pentyl gallate, present better inhibitory action in relation to antiplasmodial activity than compounds with branched structures. When comparing compounds 9 and 10, there was a significant difference (respectively 206.35 ± 46.4 µg/mL, 81.19 ± 8.54 µg/mL); and hexyl p-coumarate showed the least antiplasmodial activity against P. falciparum [60]. In compounds 11 and 12 the substituent is a substituted aromatic ring; the antiplasmodial potency of these compounds was greater than the majority of the derivatives with carbon chains; both 4-methylbenzyl p-coumarate (11) and 4-isopropylbenzyl p-coumarate (12) presented weak antiplasmodial activity with respective EC 50 of 83.98 ± 15.34 µg/mL and 93.66 ± 14.17 µg/mL. Table 2. Potential targets of compound 9 in Leishmania braziliensis identified by the homology-based target fishing approach.

L. braziliensis Target (a)
Description ID

Computational Methods
Computational target fishing approaches rely on the similarity principle: if a query compound is similar to a reference which binds a certain receptor, then the query might also bind the same receptor. Reference ligand-receptor interactions are generally extracted from large publicly available databases. However, the information in these databases is biased against human receptor-ligand associations. For this reason, the homology-based target fishing approach described in the Methods section was implemented to identify potential Leishmania braziliensis targets for compound 9. The 16 identified potential Leishmania braziliensis targets for compound 9 are provided in Table 2. Of them, seven are annotated in the UniProt database as confirmed or predicted protein kinases, suggesting a subset of possible targets enriched with this function.
Molecular docking of compound 9 (to the 16 potential targets listed in Table 2) was thus performed as described in the Methods section. It has been shown that GSK3 inhibitors are effective in binding to ATP and substrate binding sites in both Leishmania spp. and in the human homolog [61,62]. Based on these findings, we explored GSK3 for ligand binding to these two sites in Leishmania braziliensis. TOP2 is an ATP-dependent DNA topoisomerase enzyme that contains ATP and nucleotide binding sites [63,64]. For the TOP2 docking calculations, only the ATP binding cavity was explored. The results of the consensus docking protocol are summarized in Table 3. Ten (10) of the 16 modeled targets, (for compound 9), predicted more than one possible binding mode, leading to 33 ligand-receptor complexes needing further analysis Visual inspection of these 33 predicted complexes revealed meaningful binding modes for all of them. In all cases, the ligand was inside the receptor binding site, and either interaction or blocking of functionally relevant residues was noted. Table 3. Results of the molecular docking of compound 9 to its predicted potential molecular targets in Leishmania braziliensis. Considering the same consensus protocol employed for ligand pose selection, the best docking results in GSK3 targets were obtained for the HDA2 and the ATP binding sites. However, the worst docking performance for GSK3 targets was observed for the PGFS and substrate binding sites. Interestingly, the GSK3 docking results presented large differences in scores for the ATP and substrate binding sites, which could be an indicator of selectivity towards ATP binding. Overall, the docking results did not allow for confident selection of the most probable compound 9 targets. Post-processing of predicted docking complexes using Molecular Dynamics (MD) simulations along with estimations of free energies of binding can help refine target prioritization [44,65]. Calculations were thus performed (described in the Methods section) to estimate the free energies of binding for compound 9 to its predicted targets using the MD simulations data. These simulation results are demonstrated in Table 4. Table 4. Predicted free energies of binding of compound 9 to its predicted targets. The components of the free energies of binding are also included. All energies are expressed in kcal/mol. (a) van der Waals energy, (b) electrostatic energy, (c) polar solvation energy, (d) non-polar solvation energy, (e) dispersion solvation energy, (f) total gas phase free energy, (g) total solvation free energy, (h) total free energy.
The MM-PBSA calculations reveal that due to the predicted positive free energies of binding, the complexes formed by compound 9 with HDA1, PKAC3, PKAC2, and GSK3 (SBS) are unfeasible. Though negative, the ∆Gs of binding predicted for AKR, HDA2, PGFS, PGFS2, and YARS1, range from −1.08 kcal/mol to −2.01 kcal/mol. The values predict low stability for these complexes. Slightly better free energies of binding are predicted for the ALDH2, GSK3, CBPKC, CRK7, and PKC targets. Despite stable binding, (∆G TOTAL < 0) is predicted for most of the predicted compound 9 targets, as seen in the results in Table 4, where the most probable compound 9 targets are ALDH, MPK4, and TOP2. One interesting result is that of the predicted likely compound 9 targets, the binding sites of TOP2 and the kinases (MPK4, PKC, GSK3, CBPKC, and CRK7) are all ATP binding cavities. This suggests that compound 9 might interfere with ATP binding in various proteins related to diverse cellular functions. Based on the results thus obtained, we selected ALDH, MPK4, and TOP2 as targets for additional evaluation concerning compound 9 potential anti-leishmanial mechanism of action.
The interactions between compound 9 and its most probable targets were analyzed in detail along with the 200 snapshots employed for Molecular Mechanics with Poisson-Boltzmann and Surface Area solvation (MM-PBSA) calculations. The binding modes of the ligand to these receptors as predicted are shown in Figure 2. To choose a representative binding mode for each receptor, the 200 conformations of compound 9 were clustered and the centroid mode from the most populated cluster was selected.  In the case of ALDH, the ligand binds to the substrate binding channel making contact with E128, E132, Y135, Y139, F182, A185, M186, R189, R307, S471 and F477 in more than 50% of the analyzed MD snapshots. As shown in Figure 2A, the hexyl tail occupies a sub-cavity close to the NAD cofactor, and hence completely blocks access, leading to enzyme blockage. In this predicted complex, the hydroxyphenyl ring of the compound interacts with residues at the entrance of the binding pocket. Most of the interactions observed were Van der Waals and hydrophobic. One factor adding stability to this complex is the presence of hydrogen bonds between the ligand's hydroxyl moiety and E128 of ALDH. Compound 9 is also predicted to accept hydrogen bonds from the side chain of Y139 and the backbone of R189 through its carbonyl group. It is noteworthy that the hydrophobic hexyl fragment of the ligand localizes in a region containing hydrogen bond accepting capable residues such as T256, E279, T312, and the NAD cofactor. We postulate that inclusion of hydrogen bond donor substitutions in this part of compound 9 would improve the stability of this ALDH complex.
In the complex predicted for compound 9 with MPK4 ( Figure 2B), interactions with I36, V44, A57, K59, E77, M81, I90, F109, M111, Y114, T116, A120, R120, L163 and D174 occur in the majority of the MD snapshots analyzed. The interactions are mainly Van der Waals and electrostatic in nature, and as in the complex predicted with ALDH, the ligand's hexyl moiety orients toward the bottom of the cavity and the hydroxyphenyl substituent occupies the entrance of the binding site. In most of the MD snapshots studied, the hydroxyl group of the latter is predicted to achieve hydrogen bonding with either the backbone of Y114 or with side chains of R120 and D117. Another similarity between its predicted binding pose and that of ALDH is that the hydrophobic tail of compound 9 fits closely to two hydrogen bond acceptor amino acids at the receptor: E77 and D174.
The compound 9 complex predicted for the TOP2 ATP binding site ( Figure 2C) reveals a large network contacts that includes N67, A68, N71, N95, E98, I100, I116, F117, T122, S123, S124, G136, G139, F140, G141, A142 and T190. Similar to its ALDH and MPK4 targets, the compound poses in TOP2 with its hexyl part oriented towards the bottom of the binding cavity and its aromatic substituent oriented towards the cavity entrance. The phenyl substituent is also capable of forming hydrogen bonds at the drug receptor site in the TOP2 complex, either through backbones E98 and T190 or via side chains of Q79 and N95. The cavity region where the aliphatic tail of compound 9 binds contains hydrogen bonding capable residues that include N67, S123, S124 and N138. This suggests that addition of hydrogen bond donor groups to the hydrophobic tail of the compound might increase its multitarget affinity against Leishmania braziliensis through ALDH, MPK4, and the TOP2-ATP binding site.
Of the three Leishmania braziliensis targets selected as the most likely for compound 9, the aldehyde dehydrogenase function (ALDH) has been proposed as a process which (through inhibition) could be useful for development of therapeutic agents against Leishmania spp. [66]. Furthermore, there are several reports in the literature demonstrating the importance of MPK4 to Leishmania spp. in both development and infectious capability, making it an attractive for drug discovery efforts [67][68][69]. Finally, the importance of DNA topoisomerases in the life cycle of this type of parasite, and their use as drug targets has already been established in the scientific literature [70].

Conclusions
Twelve compounds were evaluated for antiparasitic potential against L. braziliensis amastigotes, against Plasmodium falciparum, and for cytotoxicity in U-937 macrophages. Eight of the derivatives were bioactive against Leishmania braziliensis, such as hexyl p-coumarate (compound 9) (4.14 ± 0.55 µg/mL; IS = 2.72) which indicated the most potent inhibitory effect. The studies performed using molecular docking suggested that compound 9 interacts with the enzymes ALDH, MPK4, and TOP2, which are essential to Leishmania braziliensis development. Some of the compounds tested were bioactive against P. falciparum yet demonstrated little antiparasitic activity. These results may be useful to other studies involving p-coumaric acid, to help obtain ever more potent and selective derivatives against these protozoan species, which continue to compromise the quality of life of various world populations.

Reagents and Chemical Characterization
The derivatives were purified by column adsorption chromatography (CC) using silica gel 60, ART 7734-MERCK, St. Louis, Missouri, EUA. Infrared spectra were performed using FTIR spectrophotometry and 1 H and 13 C NMR spectra were obtained using a Bruker-AscendTM machine, Bruker, Bremen, Germany, operating at 400 and 100 MHz. All compounds were structurally characterized as published by Lopes et al., 2019 [44].

Preparation of Compounds 9-12 by Mitsunobu Reaction
p-Coumaric acid (0.1 g; 0.61 mmol) and alcohol (0.61 mmol) were dissolved in 2.25 mL of tetrahydrofuran. The reaction mixture was stirred under magnetic stirring at 0 • C for 30 min. Subsequently, diisopropyl azodicarboxylate (0.12 mL; 0.61 mmol) and triphenylphosphine (0.16 g; 0.61 mmol) were added. The mixture was stirred at room temperature for 48-52 h and monitored with TLC [72], as published by Lopes et al. [44]. Spectroscopic data for the compounds in this study are available in Supplementary Materials.

In Vitro Cytotoxicity
The compounds were analyzed for viability of the human pro-monocytic cell line U-937 (ATCC CRL-1593.2TM) using an MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) assay based on a previously described methodology, assessing cytotoxic activity [73]. The resulting compounds were determined by measuring mitochondrial dehydrogenase activity, adding 10 mL/well of MTT solution (0.5 mg/mL) and incubation at 37 • C for 3 h. Cell viability was defined from the quantity of formazan produced and the intensity of color (absorbance) registered through optical density (O.D.) obtained at 570 nm in a (VarioskanTM Flash Multimode Reader-Thermo Scientific, Waltham, MA, USA) spectrophotometer. As a negative control, isolated cells were cultured (without the compounds) to observe viability. Doxorubicin was used as a cytotoxicity control in a similar manner. We performed two independent tests with each concentration, in triplicate [73].

In Vitro Antileishmanial Activity
The activities of the derivatives in L. braziliensis intracellular amastigotes, transfected with the green fluorescent protein gene (MHOM/CO/88/UA301-EGFP) were evaluated. Initially, human U-937 cells at a density of 3 × 10 5 cells/mL in RPMI 1640 and 0.1 µg/mL PMA (phorbol-12-myristate-13-acetate) were dispensed in a 24-well microplate and incubated at 37 • C, and 5% CO 2 . At 72 h of incubation, the cells were infected with growing (stationary phase) promastigotes; 15 parasites per cell. The plates were incubated at 34 • C and 5% CO 2 for 3 h and the cells then washed twice with phosphate buffer solution (PBS) to eliminate non-internalized parasites. Fresh RPMI-1640 was added to each well (1 mL) and the plates were incubated again to complete the infection. At 24 h of infection, the RPMI-1640 medium was replaced by a fresh culture medium containing each compound in four serial dilutions, the highest concentration being equivalent to twice the LC 50 , and the plates were incubated at 37 • C and 5% CO 2 for 72 h. The cells were then removed from the bottom plate with a trypsin/EDTA solution (250 mg), and centrifuged at 1100 rpm for 10 min at 4 • C; the supernatant was discarded and the cells were washed with 1 mL of cold PBS and centrifuged twice at 1100 rpm for 10 min at 4 • C. After this final wash, the supernatant was discarded and the cells were suspended in 500 µL of PBS.
The cells were run in a flow cytometer (Cytomics FC 500MPL), reading at 488 nm (exciting) and 525 nm (emitting) over an argon laser, and counting 10,000 events. Infected cells were determined according to green fluorescence events (parasites). Infected cells exposed to the control drug (amphotericin B) were used as antileishmanial activity controls (positive control). Infected cells incubated in the absence of any compound or drug were used as infection controls (negative control). Non-specific fluorescence was corrected by subtracting fluorescence from unstained cells. All determinations were performed in triplicate, in at least two independent experiments [74].

In Vitro Antiplasmodial Activity
Cultures of asynchronous P. falciparum 3D7 strains were adjusted for 0.5% parasitemia and 1% hematocrit in RPMI medium enriched with 3% lipid-rich bovine serum albumin-Albumax II. Afterwards, 100 µL of the parasite suspension was dispensed to each well of a 96-well cell culture plate and subsequently exposed to 100 µL of each compound at 100, 25, 6.25, or 1.56 µg/mL. The plates were incubated for 48 h at 37 • C in an atmosphere of N 2 (90%), CO 2 (5%) and O 2 (5%). After incubation, the parasites were harvested and subjected to three 20-min freeze-thaw cycles. At the same time, 100 µL of Malstat reagent (400 µL of Triton X-100 in 80 mL of deionized water, 4 g of L-lactate, 1.32 g of Tris buffer, and 0.022 g of adenine acetylpyridine dinucleotide in 200 mL of deionized water; pH 9.0), and 25 mL of NBT/PES solution (0.16 g of nitroblue-tetrazolium salt, and 0.08 g of phenazine sulfate in 100 mL of deionized water) were added to each well of an additional 96-well plate. After freezing and thawing cycles, the culture in each of the wells on the first plate was re-suspended by pipetting, and 15 µL was removed from each well and added to the corresponding well on the additional plate (containing the Malstat and NBT/PES reagents). After one hour of incubation in the dark, the color development of the lactate dehidrogenase (LDH) reaction was read on a spectrofluorometer (Varioskan, Thermo, Waltham, MA, USA) at 650 nm. The color intensity in each experimental condition was recorded in fluorescent units (F.U). Non-specific fluorescence was corrected by subtracting F.U. from the blanks. The entire experiment was performed in triplicate in at least two independent experiments. Chloroquine (QC) was used as a positive control and the culture medium was used as a negative control [75].

Statistical Analysis
Cytotoxicity was determined according to the viability and mortality percentages obtained in each experimental condition (synthesized compounds, amphotericin B, benznidazole, chloroquine, and culture medium). Viability percentages were first calculated using Equation (1).
The antiplasmodial activity of each compound was evaluated through reduction in parasite growth, calculated according to Equation (2): % parasite inhibition = 100 − [(F.U exposed parasites ÷ F.U non exposed parasites) × 100] Data for percentages of parasite reduction or inhibition were used to calculate the EC 50 using Probit linear regression modeling [76]. Activity was graded as high, moderate or low according to EC 50 values as follows: high activity when the EC 50 was <25 µg/mL; moderate activity when the EC 50 was in the range from 25 to 50 µg/mL, and low activity when the EC 50 was >50 µg/mL. The selectivity index (SI) was calculated dividing cytotoxic activity by anti-parasitic activity using the following formula: SI = CL 50 ÷ CE 50 .

Target Selection
Compound 9 (the most active) was selected for modeling studies. The potential target set for this compound was selected using a homology-based approach. Initially, potential targets of compound 9 were predicted using the similarity ensemble approach (SEA) web server [77], and the identified subset of potential targets was used as input for a BLAST search [78], for Leishmania braziliensis (taxid:5660) proteins found in the Reference Proteins Database (refseq_protein). The BLAST search was performed using its NCBI server implementation (https://blast.ncbi.nlm.nih.gov/). The Leishmania braziliensis proteins were covered (at a sequence minimum of 70%) by the BLAST sequence alignment. Those with at least 35% identity with any query protein were regarded as potential compound 9 parasite targets.

Molecular Docking
Molecular docking followed the same protocol described in our previous publications [44,65]. The initial three-dimensional conformation of compound 9 was obtained with OpenEye Omega [79]; and AM1-BCC partial atomic charges were added to it using Molcharge [80]. Since none of the predicted targets had their 3D structures solved, homology models were generated for all the studied proteins with the SWISS-MODEL server [81]. Catalytically important cofactors, including metal ions, were added to the homology models taking their templates as reference.
Gold software was used for molecular docking [82]. Receptors' binding sites were defined from the ligands present on the templates used for homology models. Any residue located at a distance lower than 6 Å from the reference ligand was selected for the active site. In the absence of a reference ligand, the binding site was defined as any amino acid within 10 Å of a point located at the center of the cavity. For primary docking, 30 different solutions were explored for compound 9 using the CHEMPLP scoring function with the Search efficiency of Gold set to 200%. The obtained 30 different ligand conformations were rescored with the GoldScore, ChemScore, and ASP scoring functions. Furthermore, during molecular docking the side chains of the residues pointing toward the binding cavity were set to flexible.
The most probable binding modes of compound 9 (for each receptor) were selected from a consensus scoring strategy. For this, the score of each ligand conformation according to scoring function j was converted to Z-scores following the transformation: where Z i,j is the Z-score of conformation i according to scoring function j,Ś j indicates the average value of scoring function j across all predicted binding modes and std S j represents the standard deviation of these values. The final aggregated score Z i of a predicted conformation was computed as the average of its Z-scores across all the CHEMPLP, GoldScore, ChemScore and ASP scoring functions. Any predicted binding mode with Z i > 1 was selected for further calculations.

Molecular Dynamics Simulations and MM-PBSA Calculations
Amber 18 [83] was employed for molecular dynamics (MD) simulations and MM-PBSA calculations. These calculations were performed using the same protocol as in our previous publication [44]. For the MD simulations, all predicted ligand-receptor complexes underwent the same modeling steps: energy minimization, heating, equilibration, production runs, and free energy of binding prediction. The Amber force field ff14SB and gaff force fields were respectively used to parameterize amino acids and non-amino acid residues. MD simulations took place in an explicit solvent for which the complexes were enclosed in truncated octahedron boxes and solvated with TIP3P water molecules. To neutralize any excess net charge, either Na + of Cl − ions were added to the solvated systems.
The prepared systems were energy minimized in two steps. First, 500 steps of the steepest descent method followed by 500 cycles of conjugate gradient were executed while constraining all atoms except solvent with a force constant of 500 kcal/mol·Å 2 . Then, all atoms were unconstrained for the 500 steps of the steepest descent algorithm followed by 1000 cycles of conjugate gradient. Long range electrostatic interactions were treated using the particle mesh Ewald (PME) method implemented in Amber, and a cutoff distance of 10 Å during all minimization steps.
The systems were then heated from 0 K to 300 K, keeping their volumes constant at 20 ps. The solute was constrained with a force constant of 10 kcal/mol·Å 2 and temperature was controlled using a Langevin thermostat with a collision frequency of 1.0 ps −1 during heating. The PME cutoff was set to 12 Å during this modeling stage. During all MD steps from this point onwards, the bonds involving hydrogen atoms were constrained and their interactions omitted by means of the SHAKE algorithm.
The heated systems were then equilibrated for 100 ps, at constant pressure (1 bar) and temperature (300 K). For equilibration, temperature was controlled during heating, and pressure was controlled with isotropic position scaling set at a relaxation time of 2 ps. The PME cutoff for equilibration was set to 12 Å. The equilibrated systems were used as inputs for 20 production runs of 2 ns length each, using the same parameters as those used during the equilibration of the systems. Atomic velocities were randomly initialized before each production run.
Compound 9 free energy of binding to its potential receptors was estimated using MM-PBSA calculations. These were performed with the MM-PBSA.py script provided with Amber 18 [84]. The MM-PBSA calculations were performed with 200 MD snapshots (one every 200 ps) selected from 20 production runs with the ionic strength set to 100 mM.

Conflicts of Interest:
The authors declare no conflict of interest.