Binding mode knowledge of new possible anti-hypertensive compounds designed in silico using Neutral Endopeptidase (NEP) as a target.

: Arterial hypertension is a health problem that affects millions of people around the world . Particularly in Chile, according to the last health survey in 2019, 28.7% of the population had this condition, and arterial hypertension complications cause one in three deaths per year. In this work, we have used molecular simulation tools to evaluate new compounds designed in silico by our group as possible anti-hypertensive agents, taking Neutral Endopeptidase (NEP) as a target, a key enzyme in the arterial hypertension regulation at the level kidney. We use docking experiments, molecular dynamics simulation, free energy decomposition calculations (by MM-PBSA method), and ligand efficiency analysis to identify the best anti-hypertensive agent pharmacokinetic and toxicological predictions (ADME-Tox). The energetic components that contribute to the complexes stability are the electrostatic and Van der Waals components; however, when the ADME-Tox properties were analyzed, we conclude that the best anti-hypertensive candidate agents are Lig783 and Lig3444, taking Neutra Endopeptidase as a target.


Introduction
Arterial hypertension is a health problem that affects millions of people around the world. The World Health Organization (WHO) has stated that more than 9.8 million people worldwide die every year from the arterial hypertension consequences (1) . Particularly in Chile, according to the last health survey in 2019, 28.7% of the population had this condition, and one in three deaths per year is caused by arterial hypertension complications (2) . For this reason, it is interesting to increase the search for more effective drugs to fight this disease.
One of the essential pharmacological targets that have mostly focused on current research is the renin-angiotensin-aldosterone system (RAAS)(3) . This system is the primary regulator of fluid and ion balance at the kidney level (4) . The RAAS comprises a series of vital proteins in the arterial hypertension regulation in the kidney. One of these enzymes is the Neutral Endopeptidase, also As shown in Figure 1, all the ligands were oriented suitably in the NEP's active center, indicating that our docking experiment acceptably reproduces the crystallographic structure obtained from the Protein Data Bank (26) . In addition to the visual analysis, we quantified the RMSD for each ligand, taking the LigHAO as a reference. As we can see in Table 1, 71.5% of the poses analyzed (250 for each compound studied) had RMSD values lower than 2 Å. The RMSD value of 2 Å was taken as a reference for a correct or incorrect (> 2 Å) docking resolution (27,28) . Only the Lig2177 compound had RMSD values greater than 2 Å.  (26,29) Base on the docking results presented in Table 1, all the complexes studied had binding energies above 7 kcal/mol. The most negative ΔGbinding was obtained in our reference ligand LigHAO (-20.86 kcal/mol). This compound was oriented in the NEP active center so that the carboxyl group formed by the oxygens O2 and O3 presented several non-covalent interactions. For example, the oxygen O2 of this functional group shows a negative charge density interacting electrostatically with the Zn 2+ of NEP (2.09 Å). Besides, oxygen 03 exhibits hydrogen bond (Hbond) interactions with the NH1 and NH2 groups of Arg717.
Furthermore, the other carboxyl group of LigHAO presents oxygen (OTX) that also interacts with the Zn 2+ of NEP, stabilizing this complex. These interactions are also found in the Neutral Endopeptidase's crystallographic structure obtained from the Protein Data Bank (PDB id: 2YB9). These results confirm that our docking experiments had excellent reproducibility ( Figure 2A).
The second most negative binding energy was found in the complex formed by Lig6199-2YB9 with a ΔGbinding value of -10.11 kcal/mol (Table 1). This ligand was oriented in the NEP active center that the hydrophobic linear hydrocarbon skeleton was inserted into a hydrophobic pocket formed by the amino acids Arg102, Phe106, Asp107, and Tyr697, stabilizing this complex ( Figure 2F).
Something similar happened with the Lig1022. This compound presented the third most negative binding energy (-9.74 kcal/mol). Like Lig6199, this ligand did not present H-bond interactions in the docking experiments because of its structural characteristic. It has a long apolar hydrophobic hydrocarbon chain, which is oriented in a hydrophobic pocket formed by the side chains of the amino acids Arg102, Phe106, Asp107, Thr708, Asp709, His711, and Arg717. This orientation confers excellent stability to the Lig1022-2YB9 complex ( Figure 2D). The Lig2177 presented the fourth most negative binding energy of all the ligands studied. The Lig2177-2YB9 complex exhibited hydrogen bond interactions with Asn704 (3.15), Tyr697 (3. 33), and with Asp709 (3.27 and 2.79). These interactions give this complex certain stability; hence its binding energy is -8.41 kcal/mol ( Figure 2C)

Molecular Dynamics Simulation.
To study the complexes' dynamic behavior, we have performed molecular dynamics simulations to know if the docking experiments' interactions are maintained during the 50 ns of simulation time. This method will also allow us to observe if these systems remain stable over time. As a stability criterion, we have quantified the RMSD, H-bond, RMSF, and Radius of gyration (Rg) parameters, which are shown below.

Root Means Squared Deviation (RMSD).
We have obtained the RMSD parameter values during the 50 ns of molecular dynamics simulation as stability proof. As shown in Figure 3, all systems remained stable over time with RMSD values lower than 1.4 Å. It should be noted that all the systems studied remained stable after 8 ns of simulation time. The Lig2177-2YB9, with an RMSD value average of 0.919±0.057 Å resulted in the most stable complex behavior. This complex presented the fourth most negative binding energy in docking experiments (out of six systems studied), which leads us to think that there are losses of interactions compared to the docking results.
The second most stable complex taking into account the RMSD values was Lig3444-2YB9 (0.923±0.068 Å). According to the molecular dynamics simulations performed, this system had similar behavior to Lig2177 in the docking results, reaffirming the previous approach about the loss of interactions in the docking experiments.
The less stable complex of all those studied (which does not mean that its behavior is unstable) was Lig6199-2YB9 with an RMSD average value of 1.10±0.105 Å, like the complex formed by our reference ligand and the NEP, which had the second-highest RMSD value of all (1.063±0.079 Å). Both systems had the two most negative energies in the docking results, indicating that the complexes' dynamic behavior differed from that found in the docking results. To find an explanation for this behavior, We will analyze the ligands' interactions designed in silico in the Neutral Endopeptidase active center at the molecular level.

Hydrogen Bond interactions (H-bond).
To explain the molecular level of the behavior found in the RMSD parameter measurements, we have quantified the hydrogen bond (H-bond) amount and stability generated during the simulation time. Besides, we can verify if the docking experiments' interactions are preserved over time.
As shown in Figure 4, the complex formed by our reference ligand and NEP (LigHAO-2YB9) was the complex with the most H-bond interactions with an average of 3.27±0.98. From 3 ns, these interactions remained stable until the end of the simulation time. To corroborate this approach, we quantify the interactions found and observe that the hydrogen bonds formed by LigHAO-O1--HN-Arg110 and LigHAO-OX--HNArg102 had 100% occupancies. This result means that the interactions were keeping by below 3A during the 50 ns of simulation time, which was our occupation's cutoff parameter. These interactions are the ones that give the most stability to this system. The second complex with the highest hydrogen bond interaction numbers was Lig6199-2YB9 (1.16±1.48). This complex's dynamic behavior was unique since the h-bonds remained stable during the first 25 ns of simulation time, reaching certain instability from nanosecond 26 to 45. From there, it stabilized again until the end of the simulation. This behavior is reflected in the sampling standard deviation, which behaved well above the sample mean. This approach is reflected in the occupation of the hydrogen bridge interactions, were in none of the cases did they exceed 50% of this parameter. The highest occupations were found in the hydrogen bonds formed by Arg114-OE2--HO-Lig6199 (34%) and Arg110-NH--O4-Lig6199 (32%). The other complexes studied had an average number of hydrogen bridge interactions below one and with occupancies below 5%, denoting certain instability in dynamic behavior.
LigHAO-2YB9 and Lig6199-2YB9 were the complexes with the most negative binding energies in the docking experiments, agrees with the results obtained in this section. However, the results obtained with the RMSD parameter do not agree with what is shown here. What is necessary to analyze other parameters extracted from the molecular dynamics simulations such as radius of gyration and RMSF, results that we will display below.

Radius of Gyration (Rg).
The results analyzed have given a certain degree of agreement between the docking experiments and the hydrogen bridges quantification obtained from the trajectories. However, there are discrepancies in the results obtained from the RMSD parameter, so we have to analyze the radius of gyration in the simulation time.
The Rg is defined as the mean square distance of the mass of a set of atoms with a common mass center (30)(31)(32) . In other words, this parameter gives us an idea of the compaction degree of the complexes studied during the simulation time. With Rg, we can compare the influence of ligands on the dynamic behavior of neutral endopeptidase. As shown in Figure 5, all the systems studied had Rg values greater than 3.5 A, higher than other complexes consulted in the literature (30)(31)(32) . However, it is necessary to explain that the complexes formed by Lig783-2YB9 and Lig3444 -2YB9 were the complexes with the lowest Rg fluctuation during the 50 ns of molecular dynamics. This result indicates that these systems were the most compact of all those studied. This behavior can be explained (at this point in the analysis) because both complexes had the lowest RMSD values (together with Lig2177-2YB9) and the lowest standard deviation values obtained from the MD simulations.
The system with the most significant fluctuation in the Rg parameter was the Lig1022-2YB9. This complex was the one with the least amount of hydrogen bond interactions (Figure 4), and it was also the complexes that had the lowest occupancy of these interactions. It should be noted that the highest H-bond interactions occupancy found in this system was formed by Asn542-NH--O-Lig1022 with 0.54%, indicating the instability of this type of interactions, which could explain the Lig1022-2YB9 low degree compaction.

Root Means Squared Fluctuation (RMSF).
The RMSF parameter will help us understand the differences in flexibility between the amino acid residues at the molecular level that make up NEP's structure when ligands designed in silico are attached to the pocket. A high value of this parameter indicates high flexibility, which could be inferred in a greater degree of freedom of movement; however, low RMSF values indicate more restricted movements during the molecular dynamics simulation (33) .
As shown in Figure 6, specific differences exist between the NEP backbone without ligands in the active center, and the complexes studied. This difference lies in the greater amino acid flexibility in the absence of ligands. However, we observed that the amino acid residues between 520 and 670 have the least freedom of movement in the backbone. Within this sequence are the amino acids that constitute the active center of NEP (Asn542, His583, Glu584, His587, and Glu646). The systems studied had lower RMSF values concerning the backbone, indicating that the ligands' binding with NEP limits the amino acids flexibility that makes up our target protein's active center. Of all the complexes, those formed by the ligands Lig1022 and Lig6199 had the highest RMSF values, indicating that these compounds in the NEP active center make the amino acids appear more flexible. This fact explains the Lig1022-2YPB complex behavior when the Rg was analyzed, which was the system that had the most fluctuation in this parameter, being the least compact complex of all those studied in this work.

Molecular Mechanics-Poison-Boltzman Surface Area methods (MM-PBSA).
To determine the energy factors that contribute to the stabilization or destabilization of the complexes studied, we have analyzed the decomposition of energy using the MM-PBSA method (34)(35)(36)(37)(38) .
The most negative binding energy was obtained by the Lig3444-2YB9 complex (ΔGbinding= -78.99± 8.67 kcal/mol) ( Table 2). This system was the second with the lowest RMSD and the second with the highest hydrogen bond number obtained from molecular dynamics simulations. This complex has been the most stable of all those studied, considering a comprehensive analysis of the results obtained so far. Our results agree with previous works obtained with this same compound but using another M4 family metalloprotein similar to Neprilysin (24) . According to ΔGbinding, other stable complexes were Lig6199-2YB9, Lig1022-2YB9, and Lig2177-2YB9, making them good candidates anti-hypertensive drugs. When analyzing the free energy decomposition presented in Table 2, we can observe that the electrostatic component, the Van der Waals interactions, and the non-polar solvation term were the stabilizing contributions of the systems studied except for the Lid1022-2YB9 complex. In this case, the electrostatic contribution was destabilizing. Considering that H-bond is a type of electrostatic interaction, it is necessary to emphasize that this complex was the one with the less hydrogen bond numbers in the molecular dynamics simulations and low stability with an occupancy of less than 0.54 % in all cases analyzed. This fact could be the explanation for the destabilizing positive value of ΔEelect.
From Table 2, we can also observe that the Van der Waals term contributes the most to the complexes stability studied, with the most negative energy of all the energy contributions calculated using the MM-PBSA method. Let's analyze the ligands' structure in this work ( Figure  8) and the docking experiments' results ( Figure 2). We can observe that the non-polar hydrocarbon skeletons present attractive hydrophobic interactions with different amino acids in the NEP's active center. This approach agrees with the non-polar solvation term results, which also contributed positively to the stability of the complexes studied.

Ligand Efficiency metrics and ADME-Tox Properties.
One of the principal objectives of this work is to analyze which of the possible antihypertensive agents is the best candidate, minimizing the risk as much as possible. According to this goal, we have performed ligand efficiency calculations and in silico predictions of the pharmacokinetic (ADME) and toxicological (Tox) properties, considering the Lipinski (40,41) , Veber (42) , and Pfizer 3/75 (43) empirical rules. Table 3 shows the behavior of these parameters for each molecule studied. There were significant differences in the Kd parameter (dissociation constant) of our reference ligand (LigHAO) for the compounds designed in silico with the lowest Kd value of all the molecules studied. The lower Kd value indicates a strong interaction between the ligand and the protein.
Of the in silico designed ligands, Lig6199 and Lig1022 had the lowest Kd values, indicating the strong interaction with Neprilysin. These results agree with the free energy calculations by the MM-PBSA method, resulting in these two ligands with the second and the third most negative binding energy of all complexes analyzed.
The ligand efficiency index (LE) is a parameter that correlates the binding energy obtained from docking experiments and the number of the atoms of the ligand, without considering the hydrogen atoms (44) . Several researchers have suggested that for a molecule to be a good drug candidate, the LE value must be greater than 0.3 kcal/mol (44)(45)(46)(47). So in this work, we take this value as a reference. Considering the results shown in Table 3, our reference ligand (LigHAO) had the highest LE value (0.6519 kcal/mol); however, of the designed ligands, only Lig3444 and Lig783 had values somewhat higher than the reference value (0.3791 and 0.3560 kcal/mol respectively). These two molecules could be good candidates for anti-hypertensive agents; however, it is necessary to analyze other parameters to reinforce or reject this hypothesis.
The binding Efficiency Index (BEI) parameter relates to the ligand'sing energy of the thP and the ligand's molecular weight (48,49) . We have taken EIB values higher than 20 and lower than 27 as a reference value base on some drugs in the market like Bortezomid (BEI = 21)(50), with EIB values in this range . As shown in Table 3, none of the molecules studied had BEI values within this range. The ligands designed in silico had values below the reference range. They could be due to these molecules' high molecular weight, which could be a negative aspect in the design of new anti-hypertensive drugs.
Another important parameter is the Lipophilic Ligand Efficiency (LLE). This aspect relates to the binding energy obtained from docking experiments with the lipophilic power (CITA). As a reference value, we consider LLE numbers between 5 and 7 units (QUOTE). It is noteworthy that the behavior of this variable was similar to BEI. None of the molecules designed in silico had LLE values in the reference range, indicating that the compounds studied have high lipophilic power, a negative behavior for a drug candidate since it can accumulate in adipose tissue and cause repeated dose toxicity.
To corroborate the above, we performed pharmacokinetic (absorption, distribution, metabolism, and elimination (ADME)) and toxicological (Tox) predictions through calculations of different parameters, shown in Table 3. As a toxicological criterion, we compared the toxicological predictions with the empirical rules of Lipinski (40,41) , Veber (42) , and Pfizer 3/75 (43) (Table 4 and 5). Table 4. Verification of the empirical rules for the prediction of Bio-availability and Toxicity of the compounds designed in silico: Red represents the violation of any of the empirical rules, and blue the No violation of any empirical rule.

Properties
LigHAO Lig783 Lig1022 Lig2177 Lig3444 Lig6199  (43) According to the results shown in Table 3 and which were integrated into Table 4, we can observe that our reference ligand meets all the criteria contemplated in the Lipinski rule. However, it does not meet the rotatable bonds criteria in Veber's rule or the TPSA (Topological Polar Surface Area) criteria from Pfizer's rule. These results indicate that LigHAO is a very flexible molecule and exceeds the polarity range, so we recommend performing experimental tests if this compound is understood as a possible anti-hypertensive compound to know if there is a toxicity mechanism.

LR VR PR LR VR PR LR VR PR LR VR PR LR VR PR LR VR PR
Of all the compounds designed in silico, Lig783 was the only one that meets all the parameters contemplated in the empirical rules of Lipinski, Veber, and Pfizer, which could be the right candidate for an anti-hypertensive agent. The other ligand to consider is Lig3444, which complies with all the Lipinski and Veber rule parameters; however, like our reference ligand, it does not comply with the TPSA parameter of the Pfizer Rule. Therefore, it is necessary to be not conclusive with this compound without doing experiments to evaluate the possible toxicity mechanism.

Computational Details
As shown previously, our group has designed a series of ligands as possible anti-hypertensive agents (18,24) , so the fundamental objective of this work is to analyze whether these designed compounds will be good candidates, taking Neutral Endopeptidase (NEP) as the target protein. To this, we have developed a sequential computational protocol, which is shown in Figure 7.
The computational protocol designed is presented in Figure 7 and will allow us to perform a comprehensive analysis of each studied molecule. The best candidates for anti-hypertensive agents will emerge from this analysis. This protocol will be explained in detail below.

Data Set.
We have designed several compounds as possible anti-hypertensive agents in previous work using the QSAR-IN and virtual screening methods (18). This result is the starting point to evaluate if any of these ligands could have this property. All molecules are represented in Figure 8, sketched using Avogadro software version 1.2.0(51) . The optimized geometry was obtained by DFT calculation at the b3lyp/ma-def2-SVP basis set implemented in Orca 4.2.1 software package (52,53) . The full optimized geometry was checked by counting their imaginary frequencies for each ligand. The fully optimized geometries of the molecules (Figure 8) were obtained using docking experiments to examine the compounds' interactions in the Neprilysin (NEP) pocket.

Docking Experiments.
The optimized geometry (from the quantum calculation of the ligands) were used for docking experiments. The compounds were prepared at pH=7.4 using Autodock Tools (54) . The NEP X-ray crystallography structure was obtained from Protein Data Bank (PDB) (29) , whose PDB id is 2YB9, resolved at 2.40 Å(39) . This protein was prepared by the addition of all hydrogen atoms at pH=7.4. The water molecules around the protein were eliminated, except those at a distance less than 5 Å from our reference ligand LigHAO. The grid box's size was 25x25x25 Å 3 around the mass centers of the LigHAO (Heteroaryl-alanine-5-phenyl oxazole). Our reference ligands discharged from x-ray crystallography structure from PDB. The grid coordinates were x = 31,959 y=-43,612 and z = 37,509 and the cluster size was 8. Neutral Endopeptidase is a metalloprotein that contains zinc (Zn 2+ ) in its active center. This fact was maintained in all docking experiments.
All docking experiments were performed under an accurate model with the flexibility of any amino acid side chain within 3 Å of the ligand. All docking was realized with the SwissDock web server (55,56) . To analyze if our docking results were correct, the LigHAO reference ligand was re-docked using the same docking protocol of the other compounds. The best docking poses were selected using binding energy (kcal/mol), FullFitness (kcal/mol) (56) , and the positional root-meansquare deviation (RMSD) (27) .
The best energetically favorable poses, the highest value of full fitness, and lowest root-meansquare deviation of each complex were selected for molecular dynamics simulations, MM-PBSA, and ligand efficiency calculations. To verify the docking results reproducibility, we calculate the root-mean-square deviation (RMSD) between the ligand designed in silico, and LigHAO as ligand reference from the Protein Data Bank crystallographic structure. These calculations were performed using the LigRMSD server 1.0 program (57) . All docking figures were built using Pymol software version 1.8(58) .

Molecular Dynamics Simulation.
We obtained the best conformational poses for each ligand-NEP complex as input for molecular dynamics simulations from docking experiments. Each complex was placed into a water box of 20x20x20 Å using the TIP3P water model (59,60) . Topologies and parameters designed in silico were obtained by the SwissParam Web Server (61) . All molecular dynamics simulations were described using CHARMM36 and CGenFF force field for the Neutral Endopeptidase and the possible anti-hypertensive compounds (62)(63)(64)(65)(66)(67) .
The ligands-NEP complexes were submitted to 50000 steps for the energy minimization using the conjugated gradient methodology, reducing any close contact. The working temperature was 298.15 K employing the weak coupling algorithm (68) . The Van der Waals cutoff was fixed to 12 Å, and we applied a backbone constraint to all complexes using the NPT ensemble. The Particle Mesh Ewald (PME) approach (69) was used for considering the long ranges of electrostatic forces. We used the velocity Verlet algorithm with a 1.0 fs time step to solve the motion equations. All the complexes were submitted to 2.0 ns of equilibration and 50 ns of molecular dynamics simulation using the NAMD 2.13 software package (70) . The trajectories analysis and scripts were perform using the VMD software version 1.9.3 (71) .

Free Energy Calculation by means MM-PBSA methods.
This work's sequential computational protocol combines docking, molecular dynamics simulation, and MM-PBSA to study the ligands designed in silico interaction with Neutral Endopeptidase. The binding free energy was calculated using g_mmpbsa package version 5.1.2 (72) , a Gromacs tool (73) to compute the ligand-NEP free binding energy. From 50 ns of molecular dynamics simulation, we extract the last 2500 frames to compute each complex's binding free energy. So the MM-PBSA method calculates the free energy decomposition into contributions. The free energy for the NEP-ligand complexes were calculated according to the following equation: In Equation 1, ΔGcomplex corresponds to the ligand-NEP complex's energy, ΔGNEP, and Δgligand as protein and ligand energy, respectively. The following equation was used to calculate the protein's free energy, ligand, and complex separately.
In Equation 2, Gx can be G complex or GNEP, or Gligand. The Ebond represents the interactions that include bond, angle, and dihedral angle, Eelect is the electrostatic energy contribution, and Evdw is a Van der Waals energy contribution. The Gpolar represents the polar free energy contribution, which was calculated using the continuum solvent Poisson-Boltzmann (PB) model included in APBS (Adaptive Poisson-Boltzmann Solver) software version 1.4.1 (74) . The non-polar free energy contribution was calculated according to the following equation:

= +
(3) In Equation 3 γ represent the coefficient related to the solvent surface tension, which, in this work, was 0.0072 kcal/mol/ Å 2 , SASA represents the solvent-accessible surface area, with an amount of 1.4 Å, and β is a fitting parameter. We also decompose the overall binding energy per residue because we need to know every amino acid's contribution. These contributions were calculated using pythons script MmPbSaStat.py (72) .

Ligand Efficiency Calculation.
Ligand efficiency metrics consist of a series of parameters that we can use to measure the relationship between the binding energy and the molecule size (47) . These parameters significantly predict how efficient a compound will be as a possible drug (44,45,47,75) . In this work, the ligand efficiency calculations were performed through several parameters shown in Table 5.
The absorption, distribution, metabolism, and excretion (ADME) properties of all ligands designed in silico were calculated from the full optimized geometry using the SwissADME web server (78) . Also, we computed other physicochemical properties. Such as molecular weight (MW), octanol/water partition coefficient (cLogP), hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), topological polar surface area (TPSA), and rotatable bond count (RB) respectively using SwissADME web server (78) . Base on the physicochemical parameters, we can predict the toxicological properties (Tox) of our ligands. They were taking into account the Lipinski (40,41) , Veber (42) , and Pfizer 3/75 toxicity empirical rules (43) (Table 6). Table 6. Empirical rules for predicting oral availability and toxicity properties of ligands studied.

Conclusions
Arterial hypertension is one of the health problems that most affect the population worldwide. Given this disease's etiology, hypertensive patients are almost forced to increase the drug dose or change it. Many researchers have been working on designing possible anti-hypertensive agents more effective and reducing as much as possible the side effects that some drugs present on the market. How we know if the compounds designed could be good anti-hypertensive agents?. To answer this question, we successfully applied a sequential computational protocol that allowed us, through a comprehensive analysis of the results, to select which compounds previously designed in silico by our group could be good anti-hypertensive agents.
The top results we have obtained were that the ligands were oriented adequately in the NEP's active center compared to our reference ligand using the docking experiments. However, we lost interactions because of hydrogen bonding from the molecular dynamics analysis, where LIg783 and Lig3444 formed the most stable complexes. This result agrees with those obtained in the calculations of free binding energy using the MM-GBSA Method in which the complex formed by Lig3444-2YB9 was with the most negative binding energy. These two ligands could be considered good candidates for anti-hypertensive agents based on the ADME-Tox predictions' favorable results. However, this result is not conclusive, first, it is necessary to perform other experimental tests that support our result.
Author Contributions: K.M.U and EL. Conceived and designed the study. J.C.C. conceived and realized the docking experiments and Molecular Dynamics Simulation, J.C.G realized the MM-PBSA calculation, E.L. performed the Ligand Efficiency Calculation, and drafted the first version of the paper, K.M.U., conceived and realized the ADME-Tox properties and contributed to the revision process of the manuscript and submit the paper. All authors reviewed the manuscript.
Funding: This work was funded by FONDECYT Iniciación grant N°11180650.

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