Identification of Potential Antiviral Inhibitors from Hydroxychloroquine and 1,2,4,5-Tetraoxanes Analogues and Investigation of the Mechanism of Action in SARS-CoV-2

This study aimed to identify potential inhibitors and investigate the mechanism of action on SARS-CoV-2 ACE2 receptors using a molecular modeling study and theoretical determination of biological activity. Hydroxychloroquine was used as a pivot structure and antimalarial analogues of 1,2,4,5 tetraoxanes were used for the construction and evaluation of pharmacophoric models. The pharmacophore-based virtual screening was performed on the Molport® database (~7.9 million compounds) and obtained 313 structures. Additionally, a pharmacokinetic study was developed, obtaining 174 structures with 99% confidence for human intestinal absorption and penetration into the blood–brain barrier (BBB); posteriorly, a study of toxicological properties was realized. Toxicological predictions showed that the selected molecules do not present a risk of hepatotoxicity, carcinogenicity, mutagenicity, and skin irritation. Only 54 structures were selected for molecular docking studies, and five structures showed binding affinity (ΔG) values satisfactory for ACE2 receptors (PDB 6M0J), in which the molecule MolPort-007-913-111 had the best ΔG value of −8.540 Kcal/mol, followed by MolPort-002-693-933 with ΔG = −8.440 Kcal/mol. Theoretical determination of biological activity was realized for 54 structures, and five molecules showed potential protease inhibitors. Additionally, we investigated the Mpro receptor (6M0K) for the five structures via molecular docking, and we confirmed the possible interaction with the target. In parallel, we selected the TopsHits 9 with antiviral potential that evaluated synthetic accessibility for future synthesis studies and in vivo and in vitro tests.


Introduction
An epidemic began in December 2019 in Wuhan, China, in which infected people suffered from pneumonia-like symptoms, which later spread throughout the world. The main cause of the infection was found to be a new virus that has structural similarities with coronaviruses related to severe acute respiratory syndrome, therefore called SARS-CoV-2 [1,2].
According to WHO data (https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports) (accessed on 28 January 2022) globally, the number of new COVID-19 cases increased in the past week (17-23 January 2022) by 5%, while the number of new deaths remained similar to those reported during the previous week. Across the six

Molecular Docking for Obtaining and Evaluating the Pose of Selected Structures and the Pharmacophoric Model
Hydroxychloroquine (HCQ) was used as a pivot molecule, according to . HCQ bind to ACE2 with KD = (7.31 ± 0.62)×10 −7 M and exhibit equivalent suppression effect for the entrance of 2019-nCoV spike pseudotyped virus into ACE2 cells [19]. In vitro studies reported that HCQ was effective against SARS-CoV-2 at a Multiplicity of Infection (MOI) of 0.01 with a 50% effective concentration (EC50) of 4.51 μM in Vero E6 cells. All MOIs (0.01, 0.02, 0.2, and 0.8) and EC50 for HCQ (4.51, 4.06, 17.31, and 12.96 μM) was satisfactory [20]. Thus, a molecule set of antimalarial analogs of 1,2,4,5 tetraoxanes and hydroxychloroquine as a pivot was preliminarily evaluated in the molecular docking study for the ACE2 target to evaluate the binding affinity and subsequent obtainment of the pharmacophoric model. The predicted inhibitory constant (pKi) was calculated using the following Equation (1) [21]:  Figure S1) as well as the pharmacophoric characteristics that were obtained via the PharmaGist web server to obtain the pharmacophoric model, in which 14 selected molecules were used as an input file and hydroxychloroquine was added as a pivot structure with an alignment score of 65.383. Subsequently, a matrix was constructed (Table 1) with the following pharmacophoric descriptors: atoms (ATM), spatial features (FEA), hydrogen bond donor (DON), hydrogen bond acceptor (ACC) and binding affinity (BA).

Molecular Docking for Obtaining and Evaluating the Pose of Selected Structures and the Pharmacophoric Model
Hydroxychloroquine (HCQ) was used as a pivot molecule, according to . HCQ bind to ACE2 with K D = (7.31 ± 0.62) × 10 −7 M and exhibit equivalent suppression effect for the entrance of 2019-nCoV spike pseudotyped virus into ACE2 cells [19]. In vitro studies reported that HCQ was effective against SARS-CoV-2 at a Multiplicity of Infection (MOI) of 0.01 with a 50% effective concentration (EC 50 ) of 4.51 µM in Vero E6 cells. All MOIs (0.01, 0.02, 0.2, and 0.8) and EC 50 for HCQ (4.51, 4.06, 17.31, and 12.96 µM) was satisfactory [20]. Thus, a molecule set of antimalarial analogs of 1,2,4,5 tetraoxanes and hydroxychloroquine as a pivot was preliminarily evaluated in the molecular docking study for the ACE2 target to evaluate the binding affinity and subsequent obtainment of the pharmacophoric model. The predicted inhibitory constant (pKi) was calculated using the following Equation (1) [21]: pKi = 10 (∆G/1.366)  Figure S1) as well as the pharmacophoric characteristics that were obtained via the PharmaGist web server to obtain the pharmacophoric model, in which 14 selected molecules were used as an input file and hydroxychloroquine was added as a pivot structure with an alignment score of 65.383. Subsequently, a matrix was constructed (Table 1) with the following pharmacophoric descriptors: atoms (ATM), spatial features (FEA), hydrogen bond donor (DON), hydrogen bond acceptor (ACC) and binding affinity (BA).
Descriptors were analyzed using the Minitab ® v. 16 software, in which the most relevant ones were used to predict the potential antiviral activity as a function of the BA value to reduce statistical inconsistencies. The ACC showed a correlation of −0.870 (strong) compared to the other descriptors, which allows us to infer that the number of hydrogen donor groups significantly interferes in the BA responses of the selected molecules. However, the contribution of each descriptor in the process of potential antiviral activity is noteworthy, as is the case of ATM with a correlation value of −0.719, FEA of The descriptors with the strongest correlations (positive and negative) were selected for evaluation by chemometric study. Dendrograms were obtained from the HCA using Minitab ® software, see Figure 2. Confirmation of the data obtained by the Pearson correlation was facilitated by generating the pharmacophoric hypotheses of the HCA in which a correlation of binding affinity (BA) as an independent variable and the structural similarity cursor was performed in the categories: higher affinity (a) and lower affinity (b), from five molecular descriptors, atoms (ATM), Spatial Features (FEA), Hydrophobic (HYD), donor (DON), and hydrogen bond acceptor (ACC).
The statistical analysis used in this study grouped structures of similar molecules into categories (Cluster). Categories are represented by a two-dimensional (2D) diagram known as a dendrogram. Molecules are represented by the branches at the bottom of the dendrogram. The similarity between clusters is given by the length of their branches so that compounds with low similarity have long branches while compounds with high similarity have short branches. The HCA method classified the molecules into two classes (high and low binding affinity) and was based on the Euclidean distance and the incremental method with full linkage [22]. HCA technique presented a similarity dendrogram in which the molecules were classified into two classes (with higher and lower binding affinity), according to their similarities, as shown in Figure 3. incremental method with full linkage [22]. HCA technique presented a similarity dendrogram in which the molecules were classified into two classes (with higher and lower binding affinity), according to their similarities, as shown in Figure 3.  The classification in the clusters considered the structural similarity in relation to the descriptors with the highest correlations. The binding affinity property stands out, which allows us to evaluate the possible interaction in the binding site, as significant values of higher BA (blue) of the structures were observed at 3-9. The cluster (blue) presents the molecules with the highest binding affinity values, consisting of molecules 3, 4, 5, 6, 7, 8, and 9. The cluster (red) classified the molecules with the lowest binding affinity value for structures 10, 11, 13, 14, 15, and 16. Pharmacophore characteristics are essential when compared to the central molecule of the process, which has three HYD groups and three ACC groups, allowing the tracking of molecules with physical and chemical characteristics closer to those of hydroxychloroquine ( incremental method with full linkage [22]. HCA technique presented a similarity dendrogram in which the molecules were classified into two classes (with higher and lower binding affinity), according to their similarities, as shown in Figure 3.  The classification in the clusters considered the structural similarity in relation to the descriptors with the highest correlations. The binding affinity property stands out, which allows us to evaluate the possible interaction in the binding site, as significant values of higher BA (blue) of the structures were observed at 3-9. The cluster (blue) presents the molecules with the highest binding affinity values, consisting of molecules 3, 4, 5, 6, 7, 8, and 9. The cluster (red) classified the molecules with the lowest binding affinity value for structures 10, 11, 13, 14, 15, and 16. Pharmacophore characteristics are essential when compared to the central molecule of the process, which has three HYD groups and three ACC groups, allowing the tracking of molecules with physical and chemical characteristics closer to those of hydroxychloroquine ( Similarity Figure 3. Dendrogram of HCA classifying structures with higher affinity (blue) and lower affinity(red).
The classification in the clusters considered the structural similarity in relation to the descriptors with the highest correlations. The binding affinity property stands out, which allows us to evaluate the possible interaction in the binding site, as significant values of higher BA (blue) of the structures were observed at 3-9. The cluster (blue) presents the molecules with the highest binding affinity values, consisting of molecules 3, 4, 5, 6, 7, 8, and 9. The cluster (red) classified the molecules with the lowest binding affinity value for structures 10, 11, 13, 14, 15, and 16. Pharmacophore characteristics are essential when compared to the central molecule of the process, which has three HYD groups and three ACC groups, allowing the tracking of molecules with physical and chemical characteristics closer to those of hydroxychloroquine (Table 2).  Tophit 313 molecules were selected using the physicochemical descriptors tracking filter on the Pharmit web server ( Table 3). The web server, in addition to the pharmacophore model, enables the special tracking filter when the stereochemical and electronic characteristics are simple and common, which allows for a more precise search with structural diversity for the construction of a small database of chemical structures. The success of molecule virtual screening shows that potential biological activity depends on the precision and specificity of the activated pharmacophore [22]. Virtual screening through a database consisting of commercial molecules from Molport and internal molecules (real/virtual files extended from the real scaffold) as an internal 3D database prepared for the virtual tracking of any model [23]. Thus, a molecular fit was applied to the molecules selected from the virtual screening based on the pharmacophore model.
The 313 molecules obtained from the rapid tracking of pharmacophoric characteristics and the application of reduction filters were subjected to the prediction of  Tophit 313 molecules were selected using the physicochemical descriptors tracking filter on the Pharmit web server ( Table 3). The web server, in addition to the pharmacophore model, enables the special tracking filter when the stereochemical and electronic characteristics are simple and common, which allows for a more precise search with structural diversity for the construction of a small database of chemical structures. The success of molecule virtual screening shows that potential biological activity depends on the precision and specificity of the activated pharmacophore [22]. Virtual screening through a database consisting of commercial molecules from Molport and internal molecules (real/virtual files extended from the real scaffold) as an internal 3D database prepared for the virtual tracking of any model [23]. Thus, a molecular fit was applied to the molecules selected from the virtual screening based on the pharmacophore model.
The 313 molecules obtained from the rapid tracking of pharmacophoric characteristics and the application of reduction filters were subjected to the prediction of pharmacokinetic properties. The plot of polar surface area and ALogP of the molecules is shown in Figure 4. nature due to the high LogP value, and of these, only one structure showed high lipophilicity and low membrane permeability due to the high LogP and molecular weight. ADME descriptors of the molecules were calculated for drug similarity studies. Intestinal absorption and blood-brain barrier penetration were predicted by developing an ADME model using the 2D PSA and AlogP98 descriptors that include 95% and 99% confidence ellipses [24,25]. Considering the established absorption, distribution, metabolism, and excretion reference parameters [26] and the pivot, the molecules selected in the previous step were evaluated within these criteria, and 182 were selected, in which they satisfy the conditions [27][28][29]. Chemical structures with less, or preferably without, violations of these rules are more likely to be administered/available orally. The results of the pharmacokinetic prediction revealed that the most active structures followed Lipinski's rule of five (R5) for oral bioavailability. The established reference ADME and pivot molecule parameters can be seen in Tables 4 and 5; see Tables S1, S2  The screening of the molecule's results in the DS-ADME model showed that of the 313 molecules subjected to pharmacokinetic prediction, only 182 have 99% confidence levels for human intestinal absorption and penetration into the blood-brain barrier (BBB). The other molecules are outside the ellipse filter of the ADME model, which indicates their lower intestinal absorption and low BBB penetration capacity. These ellipses define regions where well-absorbed molecules are expected to be found.
The compounds' good absorption or permeation through the blood-brain barrier is measured by its LogP that must be less than five [23]. Results of pharmacokinetic screening revealed that 174 molecules followed Lipinski's rule of five for oral bioavailability. Eight (8) molecule structures are out of the ellipse models because they show lipophilic nature due to the high LogP value, and of these, only one structure showed high lipophilicity and low membrane permeability due to the high LogP and molecular weight. ADME descriptors of the molecules were calculated for drug similarity studies. Intestinal absorption and blood-brain barrier penetration were predicted by developing an ADME model using the 2D PSA and AlogP98 descriptors that include 95% and 99% confidence ellipses [24,25].
Considering the established absorption, distribution, metabolism, and excretion reference parameters [26] and the pivot, the molecules selected in the previous step were evaluated within these criteria, and 182 were selected, in which they satisfy the conditions [27][28][29]. Chemical structures with less, or preferably without, violations of these rules are more likely to be administered/available orally. The results of the pharmacokinetic prediction revealed that the most active structures followed Lipinski's rule of five (R5) for oral bioavailability. The established reference ADME and pivot molecule parameters can be seen in Tables 4 and 5; see Tables S1-S3 in Supplementary Material for extended information.  BBB, blood-brain barrier (0 (Very high penetrant); 1 (High); 2 (Medium); 3 (Low); 4 (very low) [30]; Absorption, human intestinal absorption (acceptable range: range is 0-2, where 0 is a good absorption) [28]; Aqueous solubility, (acceptable range: range is 0-3, where 3 is a good solubility) [31]; Cytochrome P450 (CYP450) 2D6 inhibition (false-non-inhibitor, true-inhibitor) [28]; PPB, plasma-protein binding (false-does not bind to plasma proteins, true-binds to plasma proteins) [32,33]; Intestinal absorption (IA).
All molecules tested in the present study exhibit hydrogen bonding and hydrophobic interactions with corresponding amino acids, according to molecular docking simulations. The pivot structure did not present violation within the reference parameters (Lipinski's rule), and this same condition was observed for all molecules, which can be exemplified by the great similarity between the tested molecules, thus corroborating the studies carried out. The USFDA (Food and Drug Administration) standard toxicity risk predictor software TOPKAT (Discovery Studio, Accelrys) locates fragments within the molecule structure that indicates a potential threat to toxicity risk [30]. Toxicological predictions results for the TopHits 9 molecules can be seen in Tables 6 and 7; see Table S4 in Supplementary Material for extended information.
TOPKAT toxicity screening results for the selected compounds showed that the studied compounds do not present a risk of carcinogenicity, mutagenicity, and skin irritation, nor of skin sensitization capacity.
Similarly, the results of the USFDA rodent carcinogenicity toxicity screening, Ames mutagenicity, were negative; the test is used globally as an initial screening method to determine the mutagenic potential of new chemicals and drugs. In all parameters, ADMET and toxicological, the selected compounds (8) indicate values and characteristics superior to those of the pivot compound. Only the molecule Molport-009-499-144 showed a similar alert to hydroxychloroquine for the Ames mutagenicity test, requiring an investigation and in silico evaluation of the prediction of tolerated dose in an animal model.
Molecules Molport-005-028-274, Mol-port-009-913-111, and Molport-002-693-933 present a mild prediction for skin irritation. TOPKAT toxicity screening results for the TopHits 9 showed that the molecules studied do not present a risk of carcinogenicity, mutagenicity, and skin irritation; however, the warning can generate complications when compared to commercial compounds (Hydroxychloroquine) in development and reproduction if ingested at high doses or long-term therapeutic use in humans (see Table 7).
The carcinogenic potency data show that the molecules and pivot were within the maximum tolerated dose for rats which caused the mortality of 50% of the investigated population (TD 50 ). The Molport-005-083-430 molecule, possessed a higher value of TD 50 ; however, in a mouse model, had values within the tolerated dose (see Table 8).

Molecular Docking for ACE2 Receptor
The cryo-electron microscopy structure of the SARS-CoV-2 Spike trimer was recently reported in two independent studies. However, an inspection of the available spike protein structure revealed incomplete modeling of the RBD, particularly for the Receptor Binding Motif (RBM) that directly interacts with ACE2 [34,35]. The general structure of SARS-CoV-2 RBD together with its subunits and constituent parts, can be seen in Figure 5.
The cryo-electron microscopy structure of the SARS-CoV-2 Spike trimer was recently reported in two independent studies. However, an inspection of the available spike protein structure revealed incomplete modeling of the RBD, particularly for the Receptor Binding Motif (RBM) that directly interacts with ACE2 [34,35]. The general structure of SARS-CoV-2 RBD together with its subunits and constituent parts, can be seen in Figure  5.
SARS-CoV-2 RBD has five antiparallel β-sheet twisted strands (β1, β2, β3, β4, and β7) with short connecting helices and loops that form the nucleus [35]. Between the β4 and β7 strands in the core, there is an extended insert containing the short β5 and β6 strands, α4 and α5 helices, and loops (see Figure 5).
Given the large contact surface between Spike's RDB domain and ACE2, to carry out docking studies at this binding site, the grid configuration was centered on the Cα of the Gln493 residue located at the interface of the interaction between Spike and ACE-2, as shown in Figure 6.  SARS-CoV-2 RBD has five antiparallel β-sheet twisted strands (β1, β2, β3, β4, and β7) with short connecting helices and loops that form the nucleus [35]. Between the β4 and β7 strands in the core, there is an extended insert containing the short β5 and β6 strands, α4 and α5 helices, and loops (see Figure 5).
Given the large contact surface between Spike's RDB domain and ACE2, to carry out docking studies at this binding site, the grid configuration was centered on the Cα of the Gln493 residue located at the interface of the interaction between Spike and ACE-2, as shown in Figure 6.
This extended insert is the RBM, which contains most of the SARS-CoV-2 contact residues that bind to ACE2 [36,37]. The N-terminal peptidase domain of ACE2 has two lobes, forming the peptide substrate binding site between them.
The docking poses of all the main molecules show that they interact in a conformation that fits them into the binding pocket of the RBM. The docking poses, along with their respective interactions, are shown in Figure 7.
The generated docking poses made it possible to observe that the ligands interact with the amino acid residues of the active site of Spike RBD (PDB ID 6M0J) around the α-helix between the Tyr449-Tyr505 amino acid residues and comprised in the β-sheet between the residues of Glu35-Asp39 amino acids. In ligands, it is possible to observe hydrophobic interactions with many residues in Leu39, Tyr449, Leu452, Phe490, and Leu492; these results agree with studies in the literature [38].
In the study of molecular docking, the interactions of potential inhibitors with the amino acid residues Tyr449, Gln493, Ser494, and Tyr505 in Spike RBD are similar to those reported in the literature [39][40][41]. The best-evaluated inhibitors in terms of binding affinity were (B) MolPort-007-913-111 (−8.540 kcal/mol) and (C) MolPort-002-693-933 (−8.440 kcal/mol), in the which interactions were like those observed in the control for residues Glu35 and Ser494, contributing to the increase in binding affinity. The unusual interactions between the inhibitors were Leu39, Tyr351, Tyr 449, Phe490, Glu494, and Tyr505, and these contributions help to stabilize the active site for Spike inactivation in the RBM domain, see Figure 7.
A heatmap of the hierarchical cluster analysis of molecules can be seen in Figure 8. The analysis was performed to select the molecule with the highest representation from each group based on structural dissimilarity. It is known that similar molecules have a similar mechanism of action [42], as there is still no known drug for the treatment of SARS-CoV-2 (ACE2 target).

In Silico Determination of Biological Activity and Molecular Docking Simulations (Mpro)
In a study by Refaey et al. (2021) [43] regarding repositioning renin inhibitors as SARS-CoV-2 main protease inhibitors, five pharmacophoric characteristics were found in the pharmacophoric model, constituting two hydrogen acceptor and three hydrophobic groups, thus, these results corroborate the data obtained in this research.
Therefore, in this study, we realize the theoretical determination of biological activity for 54 structures, and only five molecules showed the potential of protease inhibitors, see Table 9  Compound MolPort-009-219-532 presented predictions to be considered both a protease and enzymatic inhibitor, see Table 9. MolPort-005 -028-274 has shown prediction protease inhibitors with a 0.134 probability to be active (Pa). In the heatmap of the molecules selected in the ACE2 target, cluster 1, the molecule that presents a chemical structure profile with greater dissimilarity from the others is MolPort-005-131-430, with a Tanimoto index (IT) of 0.23 for the 1-methoxy -4methylbenzene fragment. The interaction of the Tyr449 residue (hydrogen bonding) with the fragment stands out in the binding affinity value (molecular docking study) for the molecules of the group. In cluster 2, the molecule MolPort-002-693-933 is observed, which presents IT of 0.31 for the 1-methoxy-3methylbenzene fragment and stands out in the binding affinity value compared to the others belonging to the group, see Figure 8.

In Silico Determination of Biological Activity and Molecular Docking Simulations (Mpro)
In a study by Refaey et al. (2021) [43] regarding repositioning renin inhibitors as SARS-CoV-2 main protease inhibitors, five pharmacophoric characteristics were found in the pharmacophoric model, constituting two hydrogen acceptor and three hydrophobic groups, thus, these results corroborate the data obtained in this research.
Therefore, in this study, we realize the theoretical determination of biological activity for 54 structures, and only five molecules showed the potential of protease inhibitors, see Table 9 (  Compound MolPort-009-219-532 presented predictions to be considered both a protease and enzymatic inhibitor, see Table 9. MolPort-005 -028-274 has shown prediction protease inhibitors with a 0.134 probability to be active (Pa).
The molecular docking validation results were considered satisfactory, in which the relative positions of the crystallographic ligand and the coupled ligand were similar (Figure 9). The RMSD between the atoms of the crystallographic ligand (Mpro) and the coupled ligand was calculated to be 1.519 Å.  The molecular docking validation results were considered satisfactory, in which the relative positions of the crystallographic ligand and the coupled ligand were similar (Figure 9). The RMSD between the atoms of the crystallographic ligand (Mpro) and the coupled ligand was calculated to be 1.519 Å. According to Gowtham et al. (2008) [44] and Hevener et al. (2009) [45], the predicted binding mode using molecular docking indicates that when the RMSD is less than 2.0 Å in relation to the crystallographic pose of a respective ligand, the validation is considered satisfactory [46,47].
In control 11b (A) the interactions observed in the docking study were also similar in *Crystallographic pose (green); Docking pose (grey).  [45], the predicted binding mode using molecular docking indicates that when the RMSD is less than 2.0 Å in relation to the crystallographic pose of a respective ligand, the validation is considered satisfactory [46,47].
In control 11b (A) the interactions observed in the docking study were also similar in molecules (B) MolPort-009-219-532 and (D) MolPort-005-060-605 in relation to the active site, and others not common in molecules (C) MolPort-004-996-519, (E) MolPort-005-028-274, and (F) MolPort-009-499-144. The interactions are located around the α-helix between the amino acid residues Met49-Glu47 and in the β-sheet at residues His163-164 and Glu166, as shown in Figure 10. the amino acid residues Met49-Glu47 and in the β-sheet at residues His163-164 and Glu166, as shown in Figure 10. The selected molecules were docking with binding affinity energies to the Mpro target, and thus, a range of −8.587 to −9.012 kcal/mol was observed. The values found for the binding affinity in the docking study of all molecules are shown in the Supplementary Material, see Table S5. In the selected chemical structures, only molecule (B), MolPort-009-219-532, showed superior and/or similar results for binding affinity (−9.012 kcal/mol) compared to the values of the controls used in the study of molecular docking (11b: −8.587 kcal/mol; Lopinavir: −9.680 kcal/mol and Ritonavir: −9.594 kcal/mol). In (B), MolPort-009-219-532 interactions with hydrogen bonds and amino acid residues His163 and Gln189 were observed, suggesting the stabilization of the formed complex and the contribution of hydrophobic interactions with residues Met165-Pro168; Pi-sulfur Met49 and electrostatic interactions withthe Glu166 residue were also seen. The other molecules used in the molecular docking study, despite showing a lower binding affinity value, showed interactions similar to those observed for control groups and were not common among residues Gln189, Gln192, Asp187 and His164.
Around the α-helix (Glu47-Leu50) of the crystallographic structure, Mpro has a conformation that allows interaction with the Met49 residue, while those located in the β-sheet bind to the His163-Glu166 amino acids on the active sites of the receptor (PDB 6M0K). The major interactions involved were of the conventional hydrogen bond type observed in residues His163 and Glu166, located in the β-sheet, and Gly143 in the loop of the macromolecule. In the α-helix, Met49, Met165, and Pro168 residues appear in processes involved in hydrophobic electronic interactions, respectively. Therefore, this study showed that the results agree with studies in the literature [49,50].
In the heatmap (Figure 11 Around the α-helix (Glu47-Leu50) of the crystallographic structure, Mpro has a conformation that allows interaction with the Met49 residue, while those located in the βsheet bind to the His163-Glu166 amino acids on the active sites of the receptor (PDB 6M0K). The major interactions involved were of the conventional hydrogen bond type observed in residues His163 and Glu166, located in the β-sheet, and Gly143 in the loop of the macromolecule. In the α-helix, Met49, Met165, and Pro168 residues appear in processes involved in hydrophobic electronic interactions, respectively. Therefore, this study showed that the results agree with studies in the literature [49,50].

Synthetic Accessibility (SA) Prediction
The molecules selected in the molecular docking study were subjected to synthetic accessibility (SA) prediction and presented chemical accessibility predicted as easy, obtaining a score above 60 for both ACE2 and Mpro; that is, the molecules are easily synthesized, as shown in Table 10, see Table S6 in the Supplementary Material for extended information. Synthetic accessibility was obtained by the AMBIT web server (ambit.sourceforge.net/ reactor.html) (accessed on 31 May 2021) [51]. AMBIT calculates the complexity parameters of a molecule and issues a score ranging from 0 to 100, where 100 is the accessibility value synthetic maximum (easy synthesis) and 0 is the minimum (greater difficulty of synthesis) [52]. At the end of the virtual screening stages, nine molecules presented a better profile for SARS-CoV-2 inhibitory potential, according to Figure 12.

Prediction of Lipophilicity and Water Solubility for Promising Compounds
A parameter, commonly logP, is used to express the liposolubility of drugs, and it becomes the key point for drug planning [53]. This property affects the ability of a molecule to decompose and to decompose in non-polar environments versus aqueous environments. The nine promising compounds showed consensus logP values spanning from 3.54 to 4.26, see Table 11.
In fact, in this study, only positive logP values in the range of 1.56 to 5.71 were found. It is worth mentioning that such positive values indicate that all molecules that are highly lipophilic meet an essential criterion for a drug candidate [54].
The promising compounds showed consensus regarding logS values in the range of −4.64 to −6.26, as shown in Table 12. In this study, only negative logS values in the range −3.91 to −8.57 were found. A logS reference value for moderate solubility is between −4 and −6, −2 to −4 indicates good solubility and values greater than −6 indicate poor solubility. Solubility in water is an important requirement for any drug candidate molecule, considering its oral or parental administration, as there are many active pharmaceutical ingredients that must be administered in small volumes [55]. Therefore, we can conclude that the pivot molecule and nine promising compounds are moderately soluble in water, and the compound MolPort-005-028-274 is poorly soluble in water.
MolPort-009-499-144 76.392 Synthetic accessibility was obtained by the AMBIT web server (ambit.sourceforge.net/reactor.html) (accessed on 31 May 2021) [51]. AMBIT calculates the complexity parameters of a molecule and issues a score ranging from 0 to 100, where 100 is the accessibility value synthetic maximum (easy synthesis) and 0 is the minimum (greater difficulty of synthesis) [52]. At the end of the virtual screening stages, nine molecules presented a better profile for SARS-CoV-2 inhibitory potential, according to Figure 12.

Prediction of Lipophilicity and Water Solubility for Promising Compounds
A parameter, commonly log P, is used to express the liposolubility of drugs, and it becomes the key point for drug planning [53]. This property affects the ability of a molecule to decompose and to decompose in non-polar environments versus aqueous environments. The nine promising compounds showed consensus logP values spanning from 3.54 to 4.26, see Table 11.
In fact, in this study, only positive logP values in the range of 1.56 to 5.71 were found. It is worth mentioning that such positive values indicate that all molecules that are highly lipophilic meet an essential criterion for a drug candidate [54].     Initially, hydroxicloquine, chloroquine, and 14 structures of 1,2,4,5 tetraoxane analogues were selected from the literature with proven in vitro testing against malaria caused by Plasmodium falciparum-Sierra Leone clone D-6 to form the training set ( Figure S1, see Supplementary Materials), followed by chemometric analysis studies [56]. The molecules were optimized by the computational method DFT B3LYP 6-31G** to obtain bioactive conformation and later used as input files (in .mol and .sdf formats). Hydroxychloroquine was used as a control molecule because it has selective antimalarial activity, and the molecules were subjected to a molecular docking study in order to evaluate the binding affinity at the binding site in the receptor-binding domain in the Spike of SARS-CoV-2 linked to ACE2 with PDB ID 6M0J using the DockThor program in order to select the best bioactive pose (conformation + orientation) at the binding site for future analyses to obtain the pharmacophoric model. The methodological step can be consulted in more detail in Sections 2.2 and 3.5.

Generation and Evaluation of the Pharmacophoric Model
The input file with the pivot and molecules with the best binding affinity values were submitted to the PharmaGist web server (https://bioinfo3d.cs.tau.ac.il/PharmaGist/) (accessed on 6 March 2021) to determine the pharmacophoric characteristics: Atoms (ATM), Spatial Features (SF), Features (F), Aromatic (ARO), Hydrophobic (HYD), Acceptors (ACC) Donors (DON) [57]. An alignment score was used to choose the model and was later evaluated by the incremental method via Hierarchical Cluster Analysis (HCA) and Pearson correlation (pharmacophoric characteristics with binding affinity of studied molecules).

Selection of Molecules in the Database
The pharmacophoric model with the best alignment score was submitted to the Pharmit web server (http://pharmit.csb.pitt.edu/search.html) (accessed on 6 March 2021) for selection of the Top2000 molecules in the MolPort ® database (~7.9 million compounds) (Riga, NY, USA), based on pharmacophoric characteristics and filter (maximum and minimum) values of molecular descriptors, to increase the structural diversity in the virtual strategy [58].

Prediction of Pharmacokinetic and Toxicological Properties
Calculations of predictions of absorption, distribution, metabolism, excretion, and toxicity (ADMET) were performed using Discovery Studio v16, San Diego, CA, USA (2013) software [59,60]. These properties are important in determining the compound's success for human therapeutic use. Some important chemical descriptors correlate well with ADMET properties, such as Polar Surface Area (PSA) as a primary determinant of fraction absorption and low Molecular Weight (MW) for oral absorption. The distribution of compounds in the human body depends on factors such as the blood-brain barrier (Log BB), permeability such as Caco-2 apparent permeability, MDCK cell apparent permeability, Log Kp for skin permeability, the volume of distribution, and binding to plasma proteins (Log Khsa for protein binding).
Toxicity prediction tests were performed using Discovery Studio v.16 software via the Toxicity Prediction function by Komputer Assisted Technology (TOPKAT). Toxicity parameters included carcinogenicity in rodents, mutagenicity, the Ames test, skin irritation, eye irritation, aerobic biodegradability (AB), oral toxicity in rats (LD 50 in g/kg body weight), and whether the molecule was non-carcinogenic, non-mutagenic or non-degradable.

Molecular Docking for ACE2 Receptor with DockThor
Correct assignment of protein and ligand protonation/tautomeric states is crucial to the binding mode and its affinity predictions, requiring careful inspection of the structures. In this research, the complexes were prepared using the PDB2QR web server (https:// server.poisonboltzmann.org/pdb2pqr) (accessed on 6 March 2021) [61,62]. The assignment of protonation and tautomeric states of the ligands was performed with the Discovery Studio program, while the hydrogen atoms of the protein were added with PROPKA using pH 7.
The crystal structure of the Spike receptor binding domain of SARS-CoV-2 linked to ACE2 (Homo sapiens organism) with PDB ID 6M0J [63], resolution of 2.45 Å, and elucidated by the X-ray diffraction method.
The DockThor program uses a topology file for the ligand and cofactors (.top) and a protein-specific input file (.in) that contains the atom and partial charge types of the MMFF94 force field, both of which are generated using the built-in tools MMFFLigand and PdbThorBox. The PdbThorBox program is used to define the protein atom types and the partial charges of the MMFF94 force field. Thus, in the DockThor program, protein and ligands (including cofactor molecules) are treated with the same force field in the docking experiment [64].
The grid box configuration of each complex was automatically determined according to the reference binder when available: (1) The center of the coordinates was defined as the center of the coordinates of the ligand. (2) The grid size was defined as the largest value of the ligand axis but with a tolerance of 6 Å in each dimension. (3) Discretization (i.e., spacing between grid box points) was set to the default value of 0.25 Å.
The default parameters of the algorithm were defined as follows: (1) 24 docking runs, (2) 1,000,000 evaluations per docking run, and (3) population of 1000 individuals [65]. The quality of the protein-ligand docking score was evaluated based on the Root Mean Square Deviation (RMSD) between the best score of the docking pose and the experimental binding mode of the crystal ligand. The literature describes the common limit used to consider a highly flexible ligand coupling pose as an active type of conformation when the backbone RMSD value is ≤2.0 Å [46].

In Silico Determination of Biological Activity and Molecular Docking Simulations (Mpro)
Predictions of biological activity were performed using the online PASS web server, available at http://www.pharmaexpert.ru/passonline (accessed on 3 June 2021) [66]. Using PASS, it was possible to discover the effects of a compound based entirely on the molecular formula using MNA (multilevel neighbors of atoms) descriptors, suggesting that the biological activity is in the function of its chemical structure [67]. The drug-likeness calculations were carried out in Molinspiration analyses. Only molecules with protease inhibitors and enzyme inhibitors were selected at this stage.

Molecular Docking for Mpro Receptor
The crystal structure of the main protease (Mpro) of COVID-19 in a complex with the inhibitor 11b, PDB ID 6M0K [68] and a resolution of 1.50 Å were downloaded in the Protein Data Bank (PDB) in the format (.pdb) to perform an interaction study and receptor-ligand binding mode in the study of molecular docking. The hydroxychloroquine and 11b ligands were used as positive controls, and all water molecules and cofactors were deleted.

Structural Similarity and Synthetic Accessibility (SA) Prediction
Hierarchical clustering methods were used to select the molecules with the ChemMine Tools server, in which the measures of structural similarity of the clusters were calculated from atomic descriptors between each molecular pair, which generated a similarity matrix based on unique and common characteristics observed between molecules using the Tanimoto Index (0 = less similar and 1 = greater similarity). In the subsequent grouping steps, the similarity matrix was converted to a distance matrix by subtracting the similarity values from 1. The similarity search by ChemMine Tools allowed the structural comparison of ligands and their grouping according to similarity based on the Tanimoto Index [69].
The prediction of the synthetic accessibility (SA) of the molecules was performed using the AMBIT program (http://ambit.sourceforge.net/reactor.html) (accessed on 3 June 2021). The model for SA uses four weighted molecular descriptors, which represent different structural and topological features combined in an additive scheme [51,70]. In each target molecule or set of molecules, the algorithm calculates the molecular complexity; the stereochemical complexity is the complexity due to the presence of fused and bridged systems. The SA is issued as a score ranging from 0 to 100, where the value 100 is the maximum synthetic accessibility; that is, the molecule is more easily synthesized.

Prediction of Lipophilicity and Water Solubility for Promising Compounds
Promising molecules were evaluated with SwissADME software [53,55,56] for the prediction of lipophilicity and water solubility expressed by means of values of logP and logS, respectively. SwissADME provides five methods to predict logP values: iLOGP, xLOGP3, WLOGP, MLOGP, and Silicos-IT. iLOGP is an internal physical method of SwissADME, based on free solvation energies in 1-octanol and water calculated by the Generalized-Born model and access to surface area solvent (GB/SA). It has a performance equal to or greater than six well-established predictors.

Conclusions
The COVID-19 related pandemic is a fight that still needs to be fought by humanity, and beyond prevention by vaccination, the only way out is through the discovery of new drugs. Our study identified some potential candidates that can be used for the inhibition of Spike protein and Mpro in COVID-19.
ADMET studies have revealed that most molecules have good absorption properties and low acute toxicity values. Molecular docking studies confirmed the binding of molecules at the ACE2 active site, in which the molecule MolPort-007-913-111 had the best binding affinity value of −8.540 Kcal/mol, followed by MolPort-009-219-532 −9.012 Kcal/mol to Mpro. The similarity analysis was developed to guide future studies of molecular dynamics (MD) in the selection of these compounds. Promising compounds present binding affinity values with non-significant differences, and it may happen that when submitted to MD simulations, the behavior within the system is equivalent because they present great structural similarity. In this way, the structures of the compounds were evaluated based on the differences evident in each comparison cluster, both for the set of compounds selected in ACE2 and in Mpro. Additional experimental studies (in vitro and in vivo) need to be carried out to test possible candidates since they are easy to be synthesized, and thus better clarify the mechanism of action of the virus in the human organism.