Selene-Ethylenelacticamides and N-Aryl-Propanamides as Broad-Spectrum Leishmanicidal Agents

The World Health Organization classifies Leishmania as one of the 17 “neglected diseases” that burden tropical and sub-tropical climate regions with over half a million diagnosed cases each year. Despite this, currently available anti-leishmania drugs have high toxicity and the potential to be made obsolete by parasite drug resistance. We chose to analyze organoselenides for leishmanicidal potential given the reduced toxicity inherent to selenium and the displayed biological activity of organoselenides against Leishmania. Thus, the biological activities of 77 selenoesters and their N-aryl-propanamide derivatives were predicted using robust in silico models of Leishmania infantum, Leishmania amazonensis, Leishmania major, and Leishmania (Viannia) braziliensis. The models identified 28 compounds with >60% probability of demonstrating leishmanicidal activity against L. infantum, and likewise, 26 for L. amazonesis, 25 for L. braziliensis, and 23 for L. major. The in silico prediction of ADMET properties suggests high rates of oral absorption and good bioavailability for these compounds. In the in silico toxicity evaluation, only seven compounds showed signs of toxicity in up to one or two parameters. The methodology was corroborated with the ensuing experimental validation, which evaluated the inhibition of the Promastigote form of the Leishmania species under study. The activity of the molecules was determined by the IC50 value (µM); IC50 values < 20 µM indicated better inhibition profiles. Sixteen compounds were synthesized and tested for their activity. Eight molecules presented IC50 values < 20 µM for at least one of the Leishmania species under study, with compound NC34 presenting the strongest parasite inhibition profile. Furthermore, the methodology used was effective, as many of the compounds with the highest probability of activity were confirmed by the in vitro tests performed.


Introduction
The World Health Organization classifies Leishmania as one of the 17 "neglected diseases" that, burden primarily tropical and sub-tropical climate regions [1]. These neglected diseases are caused by parasitic agents and are considered endemic in low-income populations. To date, they have been identified in 149 countries and are responsible for anywhere between 500,000 and 1,000,000 cases annually [2].
This infectious pathology comes from the Trypanosomatidae Leishmania sp protozoa [3,4] being transmitted to humans through the bite of infected female sand flies [5]. Among the existing Leishmania species, Leishmania infantum is a flagellated protozoan that causes a border was placed, in which the pIC50 value of the withdrawn compounds comprised the range of 4.78 to 4.71, totaling 32 withdrawn compounds and 496 remaining compounds. Thus, the 248 compounds classified as active presented a pIC50 ≥ 4.79, whereas the inactive compounds had pIC50 values ≤ 4. 71.
The data set for L. infantum in its Promastigote form contained 265 compounds in total, with compounds that presented pIC50 values ≥ 4.84 (134 compounds) being active, and those with pIC50 values < 4.84 (131 compounds) being inactive.
The data set for L. braziliensis in its Promastigote form contained 265 compounds in total, with active pIC50 values ≥ 4.31 (135 compounds) and inactive pIC50 values < 4.30 (130 compounds).
The data set for L. major in its Promastigote form, with 619 compounds in total, featured 312 active compounds (pIC50 values ≥ 5.10) and 307 inactive compounds (pIC50 values < 5.08).
The Knime 3.6.2 software (Knime 3.6.2, Copyright Miner, de Konstanz Information, Zurich, Suíça, www.knime.org, accessed on 10 November 2021) was used to build and evaluate the predictive QSAR models. Due to the success of previous studies carried out in the laboratory [26,27], we chose to use 3D QSAR analysis. For this, all compounds were converted into 3D structures, saved in SDF format, and imported into Dragon 7.0 software (Kode Chemoinformatics srl, Pisa, Italy) [28] and VolSurf (Volume and Surface) [29] to obtain the descriptors. The Random Forest (RF) algorithm was used to build predictive models. The applicability domain was calculated according to the Euclidean distances present in the surveyed chemical space [27]. External cross-validation was performed to estimate the predictive power of the developed models. The performance of external models was evaluated by ROC analysis. The models were also analyzed using the MCC confusion matrix to evaluate the model globally. The complete description of the procedure performed in the Knime 3.6.2 software is shown in Figure 2.

VolSurf Descriptors
The Principal Component Analysis (PCA) was performed to identify the most significant variables based on the structures of the compounds. For the generation of the PCA models, the OH2 probe (hydrophilic probe), DRY (hydrophobic probe), N1 (amide nitrogen-hydrogen-bond donor probe) and O (carbonyl oxygen-hydrogen-bond acceptor probe) were selected [30]. The selection of the GRID force field was also performed, calculating sites of energetically favorable interactions [31]. In addition to the PCA analyses, Consensus experiments of the Principal Components Analysis (CPCA) were used to graphically evaluate the most significant blocks of descriptors, as well as the validation statistics of the variances of the first principal components (PC).

VolSurf Descriptors
The Principal Component Analysis (PCA) was performed to identify the most significant variables based on the structures of the compounds. For the generation of the PCA models, the OH 2 probe (hydrophilic probe), DRY (hydrophobic probe), N1 (amide nitrogen-hydrogen-bond donor probe) and O (carbonyl oxygen-hydrogen-bond acceptor probe) were selected [30]. The selection of the GRID force field was also performed, calculating sites of energetically favorable interactions [31]. In addition to the PCA analyses, Consensus experiments of the Principal Components Analysis (CPCA) were used to graph-ically evaluate the most significant blocks of descriptors, as well as the validation statistics of the variances of the first principal components (PC).

Homology Modelling
Sequences of enzymes and species selected in the study were obtained from the GenBank Database (National Center for Biotechnology Information, Bethesda MD, USA, https://www.ncbi.nlm.nih.gov/, accessed on 15 January 2022) [32,33], while the model structures were obtained from the Protein Data Bank library (PDB, https://www.rcsb. org/pdb/home/home.do, accessed on 15 January 2022) [34,35]. Twelve enzymes were constructed through homology models. Enzyme models were built using the molecular homology method in MODELLER 9.20 software (copyright © 2020-2020 Andrej Sali, maintained by Ben Webb in the Departments of Biopharmaceutical Sciences and Pharmaceutical Chemistry, and California Quantitative Biomedical Research Institute, Mission Bay Byers Hall, University of California, San Francisco, San Francisco, USA) [36][37][38]. Five models were generated, and the lowest energy model was chosen. The stereochemical qualities of the model were evaluated by the PSVS (Protein Structure Validation Software suite) web server (http://psvs-1_5-dev.nesg.org/, accessed on 15 January 2022), using PROCHECK [39]. PROCHECK software generates a Ramachandran graph [40], which determines which amino acid regions are allowed and not allowed in the main chain. Structural quality was assessed using VERIFY 3D software (save @ 2020 -DOE-MBI Services, http://services.mbi.ucla.edu/SAVES/, accessed on 15 January 2022) [22], and the compatibility between the protein sequence and its 3D structure, based on the chemical environment, was analyzed using WHAT IF (http://swift.cmbi.ru.nl/servers/html/index.html, accessed on 15 January 2022) [41]. Homology modeling was carried out because, with the exception of L. infantum and L. major which have two 3D structures obtained by xray crystallography, the other species of Leishmania under study (L. amazonensis, L. major and L. braziliensis) do not present all available mechanisms, so the homologous proteins were constructed.
The homology results can be found in Section S1 of Supplementary Material.

Molecular Docking Studies
Proteins were obtained from the Protein Data Bank (PDB) library (https://www.rcsb. org/, accessed on 1 November 2022) [35]. The targets were selected through bibliographical research on the mechanism of action involved in the inhibition of the parasites under study, as well as considering their structural similarity. Redocking was performed prior to validate the docking simulation. Both procedures were performed using Molegro Virtual Docker (MVD) v.6.0.1 software [42]. Enzymes and compounds were prepared according to predefined parameters in the software. In the coupling procedure (linker-enzyme), a Grid of 15 Å radius and 0.30 resolution was used, which involved the site of the binding site, defined by a known ligand for each enzyme. A model was generated to perform and evaluate the fit with expected characteristics between the ligand and the enzyme, using the MOLDOCK Score (GRID) algorithm with the scoring function and search algorithm, corresponding to Moldock. The MolDock scoring function further improves these scoring functions with a new hydrogen bonding term and new charge schemes. The docking scoring function, Escore, is defined by the following energy terms: The visualization of the established interactions was performed in the Discovery Studio Visualizer program, Biovia, 2020 (https://www.3dsbiovia.com/, accessed on 22 April 2022) [43].
Molecular Docking results can be found in Section S2 of the Supplementary material.

Predicting In Silico ADMET Properties
Calculation of ADMET parameters and toxicity risk predictions were performed using the OSIRIS Property Explorer software (Idorsia Pharmaceuticals Ltd., Allschwil, Switzerland, https://www.organic-chemistry.org/prog/peo/, accessed on 15 January 2022). The compounds were evaluated for their risk of presenting toxicity in three parameters: mutagenicity, skin irritability and reproductive toxicity. The risk assessment is carried out qualitatively so that the compounds are classified as either low risk or high risk [44,45]. For the absorption study, factors such as membrane permeability and intestinal absorption were included. We also limited our investigation to compounds that did not exceed more than two violations of Lipinski's rule and for which the logP values were in consensus.

In Vitro Assay
L. infantum, L. amazonensis, L. major and L. braziliensis were maintained in vitro as promastigotes at 26 •C in supplemented Schneiderinsectmedium (20% SFB, 100 U/mL penicillin and 100 µg/mL streptomycin, pH 7), as described by Rodrigues et al. [46]. The growth inhibition assay for the parasites was performed using promastigote forms in the logarithmic phase, which were grown in 96-well plates containing supplemented Schneider insect medium and 1 × 10 6 parasites/mL. This was performed in triplicate using different concentrations of lignans (1.56-12.5 µM), reference meglumine antimoniate drugs (200-40,000 µM), and amphotericin B (0.031-2 µM). The negative control contained neither the reference nor the tested compounds. The culture plates were maintained in a biological oxygen demand incubator (Eletrolab EL202, São Paulo, Brazil), at 26 • C, for axenic promastigotes. After 2 days under these conditions, 10 µL of 3-(4,5-Dimethylthiazol-2-yl)-2,5-Diphenyltetrazolium Bromide (MTT, 5 mg/mL) was added to each well, and the cell culture plates were incubated for 4 h before adding 50 µL of 10% sodium dodecyl sulfate (SDS) solution, to solubilize the formazan crystals. The optical density of the culture was measured in a microplate spectrophotometer reader at 540 nm (Biosystems ELx800 model, Curitiba, Brazil).

Molecular Dynamics
MD simulations were performed to estimate the flexibility of interactions between proteins and ligands using GROMACS 5.0 software (European Union Horizon 2020 Programme, Stockholm, Sweden) [47,48]. The protein and ligand topologies were also prepared using the GROMOS96 54a7 force field. The MD simulation was performed using the SPC water model of point load, extended in a cubic box [49]. The system was neutralized by the addition of ions (Cl − and Na + ) and minimized to remove bad contacts between complex molecules and the solvent. The system was also balanced at 300K, using the 100 ps V-rescale algorithm, represented by NVT (constant number of particles, volume, and temperature), followed by equilibrium at 1 atm of pressure using the Parrinello-Rahman algorithm as the NPT (constant pressure particles and temperature), up to 100 ps. DM simulations were performed in 5,000,000 steps, each at 10 ns. To determine the flexibility of the structure and whether the complex is stable close to the experimental structure, RMSD values of all Cα atoms were calculated relative to the starting structures. RMSF values were also calculated to understand the roles of nearby residues in the receptor binding site. The RMSD and RMSF graphs were generated in Grace software (Grace Development Team, http://plasma-gate.weizmann.ac.il/Grace/, accessed on 11 November 2022) and the protein and ligands were visualized in the UCSF Chimera [50]. The compound submitted to this simulation was NC34, since this molecule presented IC50 values < 10 µM in the four species of Leishmania under study, as well as the correspondence between the probability value obtained in the virtual screening and the result of the biological tests in L. infantum and L. amazonensis.

Major Signs Characterizing N-Arylpropanamide Compounds
The 16 compounds below were chosen for synthesis due to the ease of obtaining them, as well as a series of selenium compounds called selenoglycolicamides, which have similar structures have already been experimentally promising in Leishmania studies by Luis and Collaborators (2019) [16]. They were also chosen based on in silico ADME predictions in which the Selene-ethylenelacticamide derivatives showed high percentages of oral absorption. None of the compounds selected for synthesis violated Lipinski's rule, confering sufficient bioavailability.

Compounds in Study
The 77 molecules in this study were obtained through synthesis, with 48 new compounds, referred to as selene-ethylenelacticamides, derived from ethylenelactic acid . The remaining 29 compounds, N-arylpropanamides derivatives, are already reported by Sherbiny and Collaborators (2016) [20] and Huang (2021) [21] (1-29). Figure 3 displays the structure of each compound.

Quantitative Structure-Activity Relationship (QSAR) Modeling
Eight classification models were developed to perform the ligand-based virtual screening. The models used Random Forest (FR) algorithms as well as Dragon 7 [28] and VolSurf [29,30] to calculate the chemical descriptors used to analyze the activities of the compounds in L. amazonensis, L. infantum, L. braziliensis and L. major.
The RF models were evaluated for their predictive capabilities using specificity, sensitivity, positive predictive value (PPV) and negative predictive value (NPV). The performance and robustness of the models were evaluated using the Receiver Operating Characteristic Curve (ROC) and the Mathews Correlation Coefficient (MCC). Table 1 details the parameters of the four models developed with the Dragon 7 descriptors and Table 2 details the parameters found with the VolSurf descriptors.
For both of the descriptors, the cross-validation confirmed the excellent model performance, with a precision higher than 71% in all models developed. The ROC values for the depicted models were no lower than 0.8, and in consideration with the other statistics, the models were considered robust and predictive.
Given these results, the bank of selenium compounds and N-aryl-propanamide intermediates were then screened to select compounds that are predicted to be active against L. amazonensis, L. infantum, L. major and L. braziliensis. Pathogens 2023, 12, x FOR PEER REVIEW 11 of 24 (A)

Consensus Analysis
Consensus analysis was run on the generated models to increase the reliability of the independent model predictions [51,52]. After obtaining the descriptors and statistics, a consensus analysis was performed between the elaborated models. The calculation included the use of the following equation:
It was observed that both descriptors classified a few molecules similarly, and as these are analogous compounds, a pattern was elucidated. Namely, NC58 was predicted in both models above 70%, NC36 by VolSurf at 84%, and NC67 in Dragon 7 at 73%. This suggests the strong influence of the bromine atom (Br) and halogens in the structures.
The consensus of these models resulted in 38 molecules with activity above 50%. The strongest candidates within the tested applicability domain were NC58 (p = 74%), NC34 (p = 72%), NC76 (p = 71%), NC36 (p = 70%) e NC53 (p = 70%). Additionally, NC40, NC30, NC51, NC41, and NC31 showed > 61% probability upon consensus analysis. It is worth noting that molecules with electron-withdrawing groups (EWGs) in the benzene ring closest to the selenium atom (such as molecule NC58) have greater probabilities of activity so long as it did not also have an EWG in the farthest ring. We hypothesize this is due to the lower contribution of aromatic amides in the presence of EWGs.
The two descriptors showed a different classification pattern in this species of Leishmania, and it was observed that for the VolSurf descriptors the molecule with the highest probability was NC66 (p = 81%), a selenium compound with bulky (almost globular), more electronegative groups in its structure. The Dragon 7 descriptors, on the other hand, showed greater activity in the compounds of the N-aryl propanamides series, where the compound with the highest probability was NC01 (p = 82%).
In the consensus model, 55 molecules exceeded 50% probability. The strongest candidates are NC40 (p = 74%), NC53 (p = 71%), NC66 (p = 69%), NC31 (p = 68%) e NC51 (p = 68%). Other compounds of interest that were selected to move forward include NC36 (p = 68%), NC34 (p = 67%), NC19 (p = 60%), NCO2 (p = 59%), NC41 (p = 59%), and NC13 (p = 52%). Of the molecules that showed greater probability, all had a group that offered an inductive effect, except NC53, which presented a bromine atom in the ring farthest from the selenium. In this compound, the inductive effect is likely lower due to the high atomic radius of the bromine atom. Electronegative groups in the para position also reduced these inductive effects, and therefore decreased activity (NC51). Except for NC66 (whose inductive effect is decreased by two ortho electron donating groups), it also appears that activity probability increases when the benzene ring closest to the selenium was not substituted.
For this species, it was observed that the model elaborated with VolSurf descriptors was more restrictive, since it classified only 24 molecules, while the model elaborated with Dragon 7.0 descriptors classified almost twice as many compounds. Dragon 7.0 software offers a far superior number of descriptors compared to VolSurf descriptors. Despite the difference in the number of classified compounds, the importance of the methoxy groups as substituents was noticeable, since NC65 was classified as the compound with the highest probability in the two types of evaluated descriptors, corresponding to p = 77% in the VolSurf descriptors and p = 79 % in Dragon7 descriptors. From the consensus model, 25 molecules exceeded 50% activity probability. The molecules with the highest probability were: NC65 (p = 78%), NC28 (p = 73%), NC66 9 (p = 72%), NC54 (p = 65%) and NC49 (p = 63%). Compounds NC51, NC40, NC53, NC41, and NC13 delivered activity probabilities between 52% and 61% and were therefore also selected for testing.
After observing the classification patterns of the models, we noticed the importance of bulky groups in the VolSurf descriptors, since the NC50 molecule (p = 81%) presents a bis-benzyl substitution, which refers to the influence of globularity descriptors. With the Dragon 7.0 descriptors, the influence of acid groups with carboxylic acids and amine groups is noticeable, since the compound NC52 (p = 63%) was the most likely in this descriptor set, and delivers values above 60% in VolSurf descriptors.

VolSurf Descriptor Principal Component Analysis
Consensus principal component analysis (CPCA) calculates the weight of descriptor blocks and can provide insight into the most influential descriptors. The LOGS (solubility descriptors), OH2 (descriptors of hydrophilic regions), and SIZE/SHAPE (topology descriptors) were most prominent in this test set.
The PCA analysis determined that NC63, NC75, NC15, and NC27 presented structural differences compared to others in the study. The PCA corroborates the data from the CPCA, in which three groups of descriptors from the OH2 block largely explain the variability in the data through the polarity of the analyzed molecules. The most influential descriptors were identified as hydrophilic volume descriptors W1-W8, hydrophilic interaction descriptors IW1-IW4, and hydrophilic capacity descriptors CW1-CW8. Again, further analysis of these results suggests the best-performing molecules have hydrophilic regions and H-bond acceptors.

Toxicity Risk and Bioavailability In Silico Calculations
Toxicity is understood to be the inherent ability of a substance to produce harmful effects in a living organism or ecosystem. The toxic risk then, is the probability that the harmful effect, or toxic effect, occurs under the conditions of use of the substance [53,54].
Compounds with activity probabilities above 60% were evaluated for toxicity in the OSIRIS program. Parameters evaluated included mutagenicity, skin irritability and reproductive toxicity. The results are shown in Table S6, and due to their essential roles in human and animal health, must be considered in the search for new therapeutic agents [55,56]. All molecules passed these standards and were therefore maintained.
In addition to evaluating the risk of toxicity, we evaluated the rate of oral absorption and violations of Lipinski's rule, as shown in Table S7. Absorption, both transcellular and paracellular, is an important parameter when searching for drugs and medications; absorption is influenced by the molecular properties of solubility and lipophilicity [57][58][59].
Lipinski's rule, also known as the "rule of five," was one of the first parameters to correlate physicochemical properties with Absorption, Distribution, Metabolism and Excretion (ADME) parameters. This parameter presents a relationship between pharmacokinetic and physicochemical properties, indicating that a molecule will have high potential as a drug if it resembles other existing drugs. The accepted parameters include a molar mass less than 500 g/mol −1 , no more than five hydrogen bond donors or 10 hydrogen bond acceptors, and an Octanol/water partition coefficient (LogP) less than 5 [60,61].
All compounds studied had absorption rates above 60% and three compounds showed only a single violation of Lipinski's rule, these being NC70, NC73 and NC74 that present violation referring to molar mass above 500 g/mol −1 . NC69 and NC76 presented two violations of Lipinski's rule, with molecular mass greater than 500 g/mol −1 , as well as Octanol/water partition coefficient (LogP) values greater than 5, however were maintained as the pool of similar analogues displayed good absorption and bioavailability potential.

Selection of Molecules for Organic Synthesis
Following the in silico screening, 16 compounds were selected for organic synthesis and in vitro testing. The main selection criterion was a probability of activity above 50% in the consensus model, but the organic feasibility was also considered. The selected compounds were NC01, NC02, NC07, NC08, NC09, NC13, NC18, NC19, NC30, NC31, NC34, NC36, NC40, NC41, NC51 and NC53. The compounds were identified using the same ID as in silico and probabilities, and are shown in Table S8.
The demonstrated probabilities correspond to the result of the consensus of the Dragon 7.0 and VolSurf descriptors. Even though some molecules did not reach 50% activity probability in consensus, this result was maintained, as in some situations the compounds reached 50% in individual models. Similarly with the applicability domain, some molecules were maintained due to being established in the spaces comprising the domain in one of the models. Some compounds were maintained to better our understanding of the characteristics of the series under study, as well as to experimentally validate the computational methodology developed.

In Vitro Activity Assessment
The synthesized compounds were evaluated based on their ability to inhibit the promastigote form of four types of Leishmania; the results are shown in the table below (Table 3). Additionally, the results of the in vitro tests are included for comparison against the predicted percentage of activity. The performance of virtual screening studies with selenium compounds is commonly adopted in our study group, with reliable and satisfactory results, as evident in Huang's (2021) [21] selection of selenium and chlorine substituted antileishmanial compounds. In addition to the aforementioned study, the research developed by Luís et al. (2019) [16] used a combined methodology to evaluate the antileishmanial activity of selenoglycolicamide derivatives.
The results obtained from the models were corroborated by the results of the in vitro study, in which selenium compounds NC30, NC31, NC34, NC36, NC40, NC41, NC51, and NC53 presented satisfactory inhibitory concentration minimums (IC 50 ) ( Table 3). This suggests that these compounds may be able to exert leishmanicidal activity by inhibiting the promastigote form; however, further studies are needed to clearly elucidate the mechanism of action. Just the same, the in vitro results validate the developed model and recommend its reliability.

Molecular Dynamics
After analyzing the activity potential of the test compound NC34 against the four important Leishmania species under study, molecular dynamics simulations were performed with the test compound NC34 to assess the flexibility of the enzyme and the stability of interactions in the presence factors such as solvent, ions, pressure, and temperature. This information is important because it complements the results of docking and the biological tests, in addition to allowing the evaluation of the compounds that remained strongly bound to the studied enzymes in the presence of factors that are found in the host organism. For this, the following enzymes were chosen for analysis: CYP-51 from Leishmania infantum (PDB: 3L4D) (A) and Dihydroorotate Dehydrogenase from Leishmania amazonensis (B). NC34 showed a greater affinity for these proteins, corroborating the previous in silico and in vitro results of a greater inhibition for these two species. Then, the RMSD was calculated for the Cα atoms of the complexed enzyme and the structures of each ligand, separately.
For the CYP-51 target (PDB: 3L4D), the analysis of the RMSD metric of the protein ( Figure 4A) showed that the simulation stabilization for this enzyme (red line) occurred around the period of 2 ns to the period of 7ns with RMSD values corresponding to 0.27 nm. After 7ns, and until 10 ns, the enzyme presented fluctuations in the RMSD values that corresponded to 0.30 returning to 0.27 nm. With NC34 (green line), the simulation stabilization occurred around 3ns until the end of the 4ns simulation, with RMSD values corresponding to 0.35 nm and remaining stable until the total simulation time, thus demonstrating that forms greater stability when compared to the PDB linker Fluconazole. Fluconazole presented higher RMSD values compared to NC34 and in addition, it presented high peaks after the 6 ns period, demonstrating instability. The stability of the CYP-51 protein (PDB: 3L4D) is essential to keep compounds bound to the active site. For the Dihydroorotate dehydrogenase from L. amazonensis ( Figure 4B), a similar behavior was observed between the three complexes. With very close RMSD values, the three compounds began to stabilize around 1 ns. The Dihydroorotate dehydrogenase complex (red) fluctuations were observed between 5 ns and 10 ns, reaching RMSD values of 0.35 nm. For NC34, fluctuations were observed in the period of 3 to 4 ns with RMSD values of 0.34 nm and in the period of 7 ns with RMSD values of 0.37 nm. After 8 ns, there was a reduction of RMSD to 0.3 nm and reduction of fluctuations. For the control ligand (blue) fluctuations occurred at the beginning of the simulation between the period of 1 to 2 ns with RMSD values of 3.0 nm, and after this period fluctuations were observed at 5 and 9 ns with RMSD values of 0.37 nm. Thus, it can be concluded that the compound NC34 showed good stability in the target under study. The stability of the Dihydroorotate dehydrogenase protein is essential to keep the compounds bound to the active site.
When analyzing the stability of the ligands in the presence of solvents in the CYP-51 enzyme (PDB: 3L4D) ( Figure 5A), it was verified that the test compound NC34 presented higher RMSD values than the results obtained by the PDB ligand Fluconazole. Regarding the Dihydroorotate dehydrogenase ( Figure 5B), the compound NC34 showed lower RMSD values when compared to the homologous enzyme ligand. Therefore, in the presence of solvents, ions and other factors, the compound NC34 seems to be capable of establishing strong bonds with the active site of the enzyme Dihydroorotate dehydrogenase from L. amazonensis. Thus, it is suggested that NC34 tends to remain in the active site even in the presence of different factors such as temperature, pressure, solvent, and ions.
To understand the flexibility of residues and amino acids that contribute to the conformational change in CYP-51 (PDB: 3L4D) ( Figure 6A) and Dihydroorotate dehydrogenase ( Figure 6B) enzymes, root mean square fluctuations (RMSF) were calculated for each amino acid in the protein. Residues with high RMSF values suggest more flexibility and low RMSF values reflect less flexibility. Knowing that amino acids with fluctuations above 0.3 nm contribute to the flexibility of the channel structure, we identified that among the amino acids present in the CYP-51 protein (PDB: 3L4D) ( Figure 6A), the residues at positions 250, 251, 252, 253, 254 and 255 ( Figure 6B) showed RMSF values above 0.3 nm, while for the Dihydroorotate dehydrogenase protein, the residues at position 310, 311, 312 and 313 contribute to the conformational change of the protein complexed to NC34. Despite this, the residues mentioned are not components of the active site of the proteins under study, which likely contributes to NC34 remaining in the active site.

Conclusions
Sixteen selene-ethylenelactamides derivatives with the highest probability of antileishmanial activity were prioritized in silico. In addition to elucidating potential binding mechanisms, multiple physicochemical properties were calculated. The screened compounds were predicted to have high oral absorption (>60%), and were largely in accord Lipinski's rule of five, indicating good bioavailability; only seven compounds showed indications of toxicity in one or even two parameters.

Patents
The main compound of this work was patented as an unprecedented substance for tuberculosis. The process was submitted to the National Institute of Industrial Property (INPI) through the INOVA/UFPB agency, having been accepted on December 22, 2021 under case number: BR 10 2021 026090 4 and petition: 870210119497.