MOL2NET, International Conference Series on Multidisciplinary Sciences Virtual screening based on covalent docking and MM-PBSA calculations predict the drugs neratinib, sacubitril, alprostadil, trandolapril, and florbetapir as promising cruzain inhibitors useful against Chagas disease.

Graphical Abstract: Abstract. This work aimed to perform a virtual screening using a covalent docking and MM-PBSA protocol in an FDA-approved drugs library dataset to search for new compounds useful against this cruzain. Initially, 1615 FDA-approved compounds were visually inspected for the presence of chemical groups reactive against cruzain reactive cysteine (Cys 25 ), followed by the choice of the most suitable 3D structure for the virtual protocols. Thus, 241 compounds were selected and the covalent docking assays and the drugs with a fit score covalent greater than 100, were selected to the MM-PBSA calculations. Finally, the drugs neratinib, sacubitril, alprostadil, trandolapril, and florbetapir showed a covalent fit score between 102.14 and 116.59; ΔG binding values between -72.851 and -148,811 Kcal/mol calculated by MM-PBSA; and interactions with the key residues of the cruzain (Cys 25 , His 159 , Gly 23 , and Gly 65 ), showing best values than other cruzain inhibitors experimentally assayed. Our findings suggest that these drugs may be possible cruzain inhibitors, and biological assays


Introduction
Neglected tropical diseases (NTD) are groups of parasitic and bacterial disorders that affect over 149 countries of subtropical conditions [1]. These agents are more prevalent in the population living in poverty without basic sanitation, indicative of social marginalization, causing harm to countries' economies [1,2]. Public health policies related to control prevention are the main measures taken to avoid such diseases [1]. Among these, Chagas disease has been causing concern in health agencies worldwide due to the growing number of cases in recent years, endemic in more than 21 countries, with estimates of approximately 6 to 7 million infected, and 75 million in risk situation [1]. Despite presenting alarming data, there are only two drugs approved for the treatment (nifurtimox and benznidazole), being a challenge for medicinal chemists worldwide to discover a new innovative treatment [3].
Among the most explored targets in drug design in the search for new drugs against Chagas disease is cruzain [3]. This enzyme is key in several biochemical processes of the parasite, such as invasion, differentiation, proliferation, and degradation of host cell proteins, contributing to the infectious process [4]. Furthermore, cruzain is found in all stages of the parasite's life cycle (epimastigotes, trypomastigotes, metacyclic trypomastigotes, and amastigotes), becoming an excellent drug target in studies against Chagas disease [5,6].
In recent years, the strategy of discovering new drugs through drug repurposing has gained great prominence [7][8][9]. Through this method, it is possible to find new clinical indications for drugs approved for clinical use, or even drugs that already form clinical candidates, which have a clinically proven safety and efficacy profile. The advantage of repurposing a drug is that it avoids future problems related to toxicity and inadequate pharmacokinetics. This approach linked to bioinformatics, through in silico methods accelerates the discovery process and increases the probability of success in a new drug discovery campaign. [10,11].
Finally, based on the above information, a virtual drug repositioning screening campaign was carried out here using a library of 1615 FDA-approved drugs employing a covalent docking protocol and MM-PBSA calculations. From this source, it was found that the drugs neratinib, sacubitril, alprostadil, trandolapril, and florbetapir showed a covalent fit score between 102.14 and 116.59; ΔGbinding values between -72,851 and -148,811 Kcal/mol calculated by MM-PBSA; and interactions with the key residues of the cruzain (Cys 25 , His 159 , Gly 23 , and Gly 65 ), showing better values than other cruzain inhibitors experimentally assayed. Our findings suggest that these drugs may be possible cruzain inhibitors, and biological assays should be performed to confirm their potential.

Target selection and preparation
To select the best 3D structures, 18 cruzain co-crystallized with a covalent inhibitor were obtained using the database Research Collaboratory for Structural Bioinformatics Protein Data Bank database (RCSB PDB, https://www.rcsb.org/) [12]. Hydrogens were added for each of these structures, and co-crystallized ligands were removed and, subsequently, submitted to our re-docking procedure in MOL2NET, 2021, 7, ISSN: 2624-5078 https://mol2net-07.sciforum.net/ 3 a covalent approach, using GOLD ® v. 5.8.1 software [13]. The re-docking was performed by using all four algorithms, being Chemical Piecewise Linear Potential (ChemPLP), GoldScore, ChemScore, and Astex Statistical Potential (ASP). The best-fit pose was chosen for each ligand, and its value of Root-Mean-Square Deviation (RMSD) was calculated by using the PyMol ® software [14]. Finally, a heatmap was generated through Microsoft Excel ® . Thus, was chosen the most accurate structure for virtual screening simulations after our Re-docking procedures (with the lowest RMSD value). Finally, the structure 1EWO (PDB id) [15] was chosen and used in the covalent docking protocol.

Selection and preparation of ligand dataset
From the ZINC database (https://zinc.docking.org/), 1,615 FDA-approved drugs were visually inspected to the presence of chemically reactive groups against the Cys25 from cruzain. In this way, 241 were selected and, submitted to the conformational analysis using the MarvinSketch ® software [16]. Thus, ten conformations were generated for each ligand, begin choosing the conformation having the lowest energy value. The molecules were then minimized using the ArgusLab ® software [17] by applying the semi-empirical AM1 (Austin Model 1).

Covalent docking
Covalent docking studies were performed using the GOLD ® software by employing the ChemPLP algorithm (chosen after the validation procedure). Initially, was performed a conventional molecular docking, in a six Å region around the co-crystallized ligand and selected the maximum efficiency of the genetic algorithm (GA). Finally, ten binding poses were generated for each ligand, in which the most excellent fit score was chosen. With the complex obtained, the covalent bond between the electrophile group of the ligand and Cys 25 was created using Discovery Studio ® software [26]. Finally, the covalent docking using GOLD ® software was performed, selecting the electrophilic group of the ligand and the sulfur atom of Cys 25 in the same protocol described for conventional docking, generating ten complexes, in which, the most significant fit score was selected for the MM-PBSA calculations.

Molecular dynamics simulations
The trajectory files for the MM-PBSA calculations were obtained through Molecular dynamic simulations performed by using the GROMACS ® software. All complexes were obtained through molecular docking studies. Thus, water molecules were removed, while charges and hydrogens were added, using the DockPrep tool from the Chimera ® software. Then, the CHARMM36 force field was employed, followed by the TIP3P solvation method. Ligand topologies were generated by web software SwissParam (http://www.swissparam.ch/) [18]. Thus, a 1.0 nm triclinic box was created, adding water and ions at the physiological concentration (0.15M). The system was initially equilibrated in 10,000 steps by the conjugate gradient method, followed by the system's total minimization in 20,000 steps. Then, NVT (constant Number of particles, Volume, and Temperature) and NPT (constant Number of particles, Pressure, and Temperature) balances were performed at a temperature of 300 K, during ten ns. The 10 ns simulation runs were performed with the free protein and complexed with the top-hit ligands with the system assembled. Finally, the tpr and xtc files were MOL2NET, 2021, 7, ISSN: 2624-5078 https://mol2net-07.sciforum.net/ 4 used in MM-PBSA calculations. These files were used for choosing the most stable conformation, using the tool MD movie, and cluster analysis in the Chimera ® software.

MM-PBSA calculations
Molecular Mechanics/Poisson-Boltzmann Surface Area (MM-PBSA) method calculates the energy of interactions between ligand and a macromolecule. These calculations are generally applied in high-throughput virtual screenings to reduce false positives [19]. Thus, this method calculates the Gibbs free-biding energy (ΔGbinding) based on van der Waals and electrostatic (unbound) interactions between the ligand and its receptor during a molecular dynamics simulation [20]. Thus, the ΔGbinding is calculated by the difference between the free energy of complex protein-ligand (Gcpx) and unbound protein and ligand (Grec). Finally, the ΔGbinding was considerate such as summation of the changes in the solvation entropy (−TΔSsol), binding energy (ΔEbind) and conformational entropy (−TΔSconf), as shown in equation 1 [21]. MM-PBSA calculations were performed using the g_mmpbsa tool [19] through the trajectory files obtained after the molecular dynamics simulation at 100 ns, using the GROMACS ® software. Then, ΔGbinding values were determined as the average of free-interaction and solvation energies [20].

Interaction analyses of most stable conformation.
After the protocol of molecular dynamics, the most stable conformation of each complex was analyzed for their interactions with the key residues (Cys 25 , His 159 , Gly 23 , and Gly 65 ) using the Discovery Studio ® software [22]. The 3D figure interactions were generated using the Chimera ® software [23].

Screening protocol validation
For validation of screening protocol, was performed a re-docking in a covalent docking protocol of 18 crystal structures of the cruzain available in the PDB database using the four algorithms (ChemPLP, GoldScore, ChemScore, and ASP) of the GOLD® software. This procedure is according to other previously published by our research group [24][25][26]. Thus, the RMSD of the ligands was calculated and the PDB 1EWO in the algorithm ChemPLP was chosen, showing the best RMDS value (0.275) between the ligand docked and co-crystalized. Next, this crystal structure was used in a covalent docking protocol of the 18 ligands co-crystallized in available in the PDB database. The ligands showed a fit score between 49.85 -94.18 (Table 1), and a covalent fit score varying between 49.91 -141.52 (Table 1). In addition, it was observed that the compounds with conventional fit score and covalent fit score better than 55 and 90, respectively, showed interactions with the main residues related to the cruzain inhibition (Cys 25 , His 159 , Glu 23 , Gly 65 , and Glu 205 ) [3,5,6,26,27] (Table 1). Other works used a similar protocol to search new hits to drug design [10,28,29]. In this way, these values were used as a start point to filter the drugs in the virtual screening covalent docking-based (next section).

Virtual screening covalent docking-based
With PDB crystal structure, GOLD ® algorithm selected and the start point values for screening the drugs, the next stage was built the ligand dataset. For this, 1615 FDA-approved drugs from the ZINC database were visually analyzed for choosing the drugs with chemical groups that can interact covalently with the sulfur atom of Cys 25 from cruzain [3]. This analysis resulted in the selection of 241 drugs that contained ester, carbon α-halogenated, nitrile, amide, and alkenes, and alkyne. The covalent docking was performed in this dataset, resulting in the selection of 10 compounds with conventional fit score and covalent fit score better than 55 and 90, respectively ( Table 2) and interactions with the key residues ( Table 2), similar to the know inhibitors showed in Table 1.
Interestingly, our screening protocol was able to identify drugs of the most diverse therapeutical indication and with previously known anti-Trypanosoma activity. The drugs ZINC1853205 (trandolapril) and ZINC3792417 (sacubitril) are inhibitors of the angiotensinconverting enzyme (ACE), and the ZINC29319828 (eprosartan) is an angiotensin receptor blocker, in which several pieces of evidence indicate that this therapeutics class shows benefices against chronic infection of Trypanosoma cruzi [30][31][32][33][34]. Furthermore, ZINC1536109 (pralatrexate) and ZINC3916214 (neratinib) are anticancer drugs, in which evidence point that the first is a solute carrier proteins (SLCs) inhibitor, similar to pentamidine, a drug used against African trypanosomiasis, and can be useful against the Chagas diseases [35,36]. On other hand, the ZINC3916214 (neratinib) is used against HER2* breast cancer with an electrophilic group that some authors associate as cysteine proteases inhibitor [36]. In addition, the drug ZINC4098633 (polydatin) is a resveratrol glucoside, and the resveratrol is a known compound with anti-Trypanosoma cruzi activity [37,38]. Other drugs with known activity are ZINC1530575 (capsaicin), previously identified as a Trypanosoma cruzi arginine kinase (TcAK) inhibitor, and can be promising results in biological assays against this disease [39]. Finally, the drugs ZINC3813078 (epoprostenol), ZINC3813088 (alprostadil), ZINC43207237 (florbetapir) are prostacyclin analog, a prostaglandin analog, and positron emission tomography (PET), respectively, that can be explored their potential against cruzain and Chagas diseases [40][41][42].

MM-PBSA calculations
The MM-PBSA calculations were performed to re-score the best drugs identified here, and the results are shown in Table 3. The drug neratinib shows de best free binding energy, with a high ΔGbinding of -148.811 ± 15.601 kJ/mol, followed by the drugs trandolapril, florbetapir, sacubitril, and alprostadil (ΔGbinding value of -94.222 ± 15.615, -93.683 ± 15.577, -77.117 ± 15.489, and -72.851 ± 32.564 kJ/mol, respectively). The high free binding energy values indicate that the drugs can be useful against cruzain. These results are according to other works, such as Silva et al., [26], in which their most potent cruzain inhibitor shows an ΔGbinding value of -50.71 kJ/mol. In addition, the other drugs identified pralatrexate and polydatin showed great free binding energy values (ΔGbinding value of -67.445 ± 26.952, and -54.591 ± 14.666 kJ/mol, respectively), and can be evaluated to determine their potential. The drugs epoprostenol, capsaicin, and eprosartan (ΔGbinding value of -34.764 ± 43.137, -19.080 ± 58.756, and 446.749 ± 25.311 kJ/mol, respectively) showed the low values. Our findings indicate that the drugs can be investigated in biological assay and offer promising results.
In addition to the free binding energy, the drugs were analyzed concerning the main interactions involved in the formation of the drug-cruzain complex. Interesting, the drugs with the most affinity to the cruzain (neratinib, trandolapril, and florbetapir) showed the best SASA energy value (-21.582 ± 1.225, -17.836 ± 1.200, -17.967 ± 2.242 kJ/mol, respectively), indicating bigger permanence in a hydrophobic environment, away from the external environment exposed to solvent, which can be related to the most excellent affinity to the target [43]. On the other hand, the drugs epoprostenol, capsaicin showed high SASA energy values (-6.101 ± 6.184, and -4.324 ± 5.591 kJ/mol), justifying the low affinity. Also, all drugs showed the values of Van der Waals energy better than electrostatic energy (see Table 3), which indicates that van der Walls interactions are the most important to the formation of the drug-enzyme complex [44]. Finally, the polar solvation energy with positive values (see Table 3) and high electrostatic energy suggest that the polar contribution is unfavorable to the binding processes of these drugs with the cruzain [45]. Thus, these finding indicates new insight about chemical compounds that can be affinity to this target.
The top 5 hit compounds (neratinib, sacubitril, alprostadil, trandolapril, and florbetapir) were calculated the contribution energy per-residue by MM-PBSA to determine the mains related to the interaction process, shown in Figure 1. In fact, all top show the Cys 25 as the best energy (between -1.14 and -3.29 kJ/mol), suggesting the strong interactions with this residue. Besides, other residues such as Leu 67 (between -2.12 and -7.34 kJ/mol), Asp 135 (between -2.89 and -3.79 kJ/mol) showed their importance in forming the complex. Authors such as Wiggers et al. [46] highlight the importance of the interaction with these residues to the cruzain inhibition. Finally, the drug of the most affinity energy (neratinib) showed great interaction with Gln 23 (-3.58), Cys 25 (-3.24), His 159 (-3.80), Glu 205 (-5.13), and other residues, which contribute to the most excellent affinity.

Interaction analysis of top 5 hits at the active site
Using the trajectory of the molecular dynamics simulations, the cluster analyses were performed for choosing the most stable conformation, following the interaction analysis at the active site of cruzain shown in Figure 2. The drug neratinib showed strong interactions with Cys 25 (H-bond) and His 159 (C-H-bond) related to the catalytic activity of this target, justifying the most excellent affinity. In addition, the compounds interacted with other residues, such Leu 67 (VDW), Asp 158 (VDW), and Ala 133 (VDW) similar to known inhibitors [47,48]. All drugs showed interactions with Ala133, Asp158, Leu67, similar to the findings of Cardoso et al. [48], in which their compounds exhibited IC50 between 0.01 and 9.3 μM against cruzain. These results suggested that the compounds can be promising in biological assays.

Conclusions
The discovery of a new drug against the Chagas disease is one of the main efforts of medicinal chemistry. Despite the existence of treatments against this disease, the drugs showed high toxicity and in the chronic stage of infections are ineffective. Here, a virtual drug repurposing screening was able to identify new drugs that can be overcome several challenges in new treatments for Chagas disease. In fact, the drugs neratinib, sacubitril, alprostadil, trandolapril, and florbetapir showed the best results employing a protocol based on covalent docking and MM-PBSA calculations, similar to the experimental data of known inhibitors. These findings suggest new drugs that can be evaluated in vitro and in vivo assay, to prove their potential anti-Chagas disease.