Design and Identification of Inhibitors for the Spike-ACE2 Target of SARS-CoV-2

When an epidemic started in the Chinese city of Wuhan in December 2019, coronavirus was identified as the cause. Infection by the virus occurs through the interaction of viral S protein with the hosts’ angiotensin-converting enzyme 2. By leveraging resources such as the DrugBank database and bioinformatics techniques, ligands with potential activity against the SARS-CoV-2 spike protein were designed and identified in this investigation. The FTMap server and the Molegro software were used to determine the active site of the Spike-ACE2 protein’s crystal structure. Virtual screening was performed using a pharmacophore model obtained from antiparasitic drugs, obtaining 2000 molecules from molport®. The ADME/Tox profiles were used to identify the most promising compounds with desirable drug characteristics. The binding affinity investigation was then conducted with selected candidates. A molecular docking study showed five structures with better binding affinity than hydroxychloroquine. Ligand_003 showed a binding affinity of −8.645 kcal·mol−1, which was considered an optimal value for the study. The values presented by ligand_033, ligand_013, ligand_044, and ligand_080 meet the profile of novel drugs. To choose compounds with favorable potential for synthesis, synthetic accessibility studies and similarity analyses were carried out. Molecular dynamics and theoretical IC50 values (ranging from 0.459 to 2.371 µM) demonstrate that these candidates are promising for further tests. Chemical descriptors showed that the candidates had strong molecule stability. Theoretical analyses here show that these molecules have potential as SARS-CoV-2 antivirals and therefore warrant further investigation.


Introduction
The coronavirus, a non-segmented virus that acts on single strands of RNA, is a member of the family Coronaviridae and the order Nidovirales. It possesses the type III transmembrane glycoprotein and the envelope protein, which together make up the S proteins in its membrane [1]. Thirty species have previously been identified as being involved in the infection of many different animals, including humans. Many diseases caused by members of this viral family in people are mild, but two beta-type coronaviruses have worsened and spread epidemics. The first of these was SARS-CoV, which was first discovered in 2003 and led to an outbreak that killed 436 people across six countries and

Determination of Theoretical Activity Site
For this study, the protein crystal structure of the SARS-CoV-2 spike receptor binding domain bound with ACE2 (Spike-ACE2) (PDB ID: 6M0J) was selected (a design that represents the binding of the spike protein with the host biological receptor). Finding novel methods to identify the activity region was important because the target's active site has not yet been identified and lacks the complexed ligands required for methods such as redocking.
At first, the FTMap server pointed out three probable regions of hot spots that can be investigated as the potential active site of the Spike-ACE2 protein ( Figure 2). Region 2 presented a more significant number of probe molecules; however, for further confirmation, hydroxychloroquine was docked in all regions and the one in which the binding affinity was the best selected.
Hot spots R2 and R3 were predominant in the region of the ACE2 protein, whereas R1 is concentrated between the protein bindings of the spike protein with ACE2. It was hypothesized that the active region is more likely to be clustered between the two proteins, preventing selectivity for ACE2 alone. Molecular docking supported this hypothesis: Hydroxychloroquine was shown to bind to R1 more readily than in other regions. Therefore, the coordinates of the active theoretical site for this target were determined to be X = −32.770, Y = 24.940, and Z = 6.690.
To reinforce the hypothesis of the active theoretical site, the Molegro software was used to analyze the cavities present in the protein, and the presence of a cavity in the same R1 region found in the FTMap was noticed, with a difference of ≈1 Å in its coordinates

Determination of Theoretical Activity Site
For this study, the protein crystal structure of the SARS-CoV-2 spike receptor binding domain bound with ACE2 (Spike-ACE2) (PDB ID: 6M0J) was selected (a design that represents the binding of the spike protein with the host biological receptor). Finding novel methods to identify the activity region was important because the target's active site has not yet been identified and lacks the complexed ligands required for methods such as redocking.
At first, the FTMap server pointed out three probable regions of hot spots that can be investigated as the potential active site of the Spike-ACE2 protein ( Figure 2). Region 2 presented a more significant number of probe molecules; however, for further confirmation, hydroxychloroquine was docked in all regions and the one in which the binding affinity was the best selected.
Hot spots R2 and R3 were predominant in the region of the ACE2 protein, whereas R1 is concentrated between the protein bindings of the spike protein with ACE2. It was hypothesized that the active region is more likely to be clustered between the two proteins, preventing selectivity for ACE2 alone. Molecular docking supported this hypothesis: Hydroxychloroquine was shown to bind to R1 more readily than in other regions. Therefore, the coordinates of the active theoretical site for this target were determined to be X = −32.770, Y = 24.940, and Z = 6.690.
To reinforce the hypothesis of the active theoretical site, the Molegro software was used to analyze the cavities present in the protein, and the presence of a cavity in the same R1 region found in the FTMap was noticed, with a difference of ≈1 Å in its coordinates

Selection and Optimization of Structures
Based on research demonstrating the inhibitory effects of hydroxychloroquine and chloroquine on SARS-CoV-2 targets, we postulated that compounds with similar structural features and chemical properties would also be potential inhibitors for this target. In this manner, antiparasitic drugs that have already received approval from the relevant agencies were sought to act as the basis for this study.
One hundred antiparasitic drugs were chosen from the DrugBank server and put through a Tanimoto similarity search, which compares the properties of these antiparasitic medications to those of hydroxychloroquine and chloroquine. This is one of the best methods for identifying how similar certain types of molecules are to one another. Only nine of these had an index above 0.32 (Figure 4), which was sufficient for moving on to creating pharmacophoric models and carrying out the virtual screening.
Chloroquine had a Tanimoto index of 0.96 compared to hydroxychloroquine's structure, while quinacrine, piperaquine, amodiaquine, primaquine, pyronaridine, mefloquine, quinidine, and quinine had Tanimoto   typically located in the most excellent cavity pocket present, and this relationship between cavity pockets and their presence has already been examined.

Selection and Optimization of Structures
Based on research demonstrating the inhibitory effects of hydroxychloroquine and chloroquine on SARS-CoV-2 targets, we postulated that compounds with similar structural features and chemical properties would also be potential inhibitors for this target. In this manner, antiparasitic drugs that have already received approval from the relevant agencies were sought to act as the basis for this study.
One hundred antiparasitic drugs were chosen from the DrugBank server and put through a Tanimoto similarity search, which compares the properties of these antiparasitic medications to those of hydroxychloroquine and chloroquine. This is one of the best methods for identifying how similar certain types of molecules are to one another. Only nine of these had an index above 0.32 (Figure 4), which was sufficient for moving on to creating pharmacophoric models and carrying out the virtual screening.

Selection and Optimization of Structures
Based on research demonstrating the inhibitory effects of hydroxychloroquine and chloroquine on SARS-CoV-2 targets, we postulated that compounds with similar structural features and chemical properties would also be potential inhibitors for this target. In this manner, antiparasitic drugs that have already received approval from the relevant agencies were sought to act as the basis for this study.
One hundred antiparasitic drugs were chosen from the DrugBank server and put through a Tanimoto similarity search, which compares the properties of these antiparasitic medications to those of hydroxychloroquine and chloroquine. This is one of the best methods for identifying how similar certain types of molecules are to one another. Only nine of these had an index above 0.32 (Figure 4), which was sufficient for moving on to creating pharmacophoric models and carrying out the virtual screening.
Following that, these antiparasitics were put through structure optimization along with hydroxychloroquine and chloroquine to achieve their optimal conformations and prevent errors in the in silico models. Based on the technique of Ferreira et al. (2019) [14], molecular mechanics using MM+ force field through the HyperChem program was used.

Determination of Pharmacophoric Regions
Following the Tanimoto analysis, the molecules were aligned using the Discovery Studio program and transmitted to the PharmaGist web server, with the structure of hydroxychloroquine being used as a pivot. There were three pharmacophoric features produced by alignment. The coordinates of each pharmacophoric property, two aromatic (AR) and one hydrogen bond acceptor (HBC), are shown in Table 1. The search resulted in a score of 28.062 from the PharmaGist server, the molecules were aligned, taking into account their Tanimoto values, where the reference has a value of 1. As a result, a matrix of pharmacophoric characteristics, atoms (ATM), characteristics (F), spatial characteristics (SF), aromatics (ARO), hydrophobic (HYD), donors (DONN), acceptors (ACC), resulting from the alignment of molecules were generated as shown in Table 2. Chloroquine had a Tanimoto index of 0.96 compared to hydroxychloroquine's structure, while quinacrine, piperaquine, amodiaquine, primaquine, pyronaridine, mefloquine, quinidine, and quinine had Tanimoto indices of 0.78, 0.74, 0.52, 0.47, 0.42, 0.36, 0.35, and 0.35, respectively. Tafenoquine had the lowest Tanimoto indices of 0.33.
Following that, these antiparasitics were put through structure optimization along with hydroxychloroquine and chloroquine to achieve their optimal conformations and prevent errors in the in silico models. Based on the technique of Ferreira et al. (2019) [14], molecular mechanics using MM+ force field through the HyperChem program was used.

Determination of Pharmacophoric Regions
Following the Tanimoto analysis, the molecules were aligned using the Discovery Studio program and transmitted to the PharmaGist web server, with the structure of hydroxychloroquine being used as a pivot. There were three pharmacophoric features produced by alignment. The coordinates of each pharmacophoric property, two aromatic (AR) and one hydrogen bond acceptor (HBC), are shown in Table 1. Following that, these antiparasitics were put through structure optimization along with hydroxychloroquine and chloroquine to achieve their optimal conformations and prevent errors in the in silico models. Based on the technique of Ferreira et al. (2019) [14], molecular mechanics using MM+ force field through the HyperChem program was used.

Determination of Pharmacophoric Regions
Following the Tanimoto analysis, the molecules were aligned using the Discovery Studio program and transmitted to the PharmaGist web server, with the structure of hydroxychloroquine being used as a pivot. There were three pharmacophoric features produced by alignment. The coordinates of each pharmacophoric property, two aromatic (AR) and one hydrogen bond acceptor (HBC), are shown in Table 1. The search resulted in a score of 28.062 from the PharmaGist server, the molecules were aligned, taking into account their Tanimoto values, where the reference has a value of 1. As a result, a matrix of pharmacophoric characteristics, atoms (ATM), characteristics (F), spatial characteristics (SF), aromatics (ARO), hydrophobic (HYD), donors (DONN), acceptors (ACC), resulting from the alignment of molecules were generated as shown in Table 2.  The search resulted in a score of 28.062 from the PharmaGist server, the molecules were aligned, taking into account their Tanimoto values, where the reference has a value of 1. As a result, a matrix of pharmacophoric characteristics, atoms (ATM), characteristics (F), spatial characteristics (SF), aromatics (ARO), hydrophobic (HYD), donors (DONN), acceptors (ACC), resulting from the alignment of molecules were generated as shown in Table 2.

Hierarchical Cluster Analysis (HCA) and Molecular Overlay
A structure-activity relationship analysis was carried out utilizing Pearson's correlation using the physicochemical properties received from PharmaGist based on the Tanimoto index values. Among the observed characteristics, the hydrogen acceptors (ACC) presented the best correlation with the Tanimoto index with a value of −0.564, thus being an important parameter of molecular similarity, followed by the characteristics (F) that present a value of −0.547. Spatial characteristics (SF) showed −0.406, hydrophobic and donors had low values of −0.295 and −0.288, respectively (Table 3). Hierarchical cluster analysis (HCA) was performed with the help of Minitab 15 software ( State College, PA, USA). A dendrogram of the pharmacophoric characteristics was generated ( Figure 5), which confirmed the Pearson correlation values, demonstrating a greater approximation of ACC with the Tanimoto index (TI) and a departure from the HYD and DONN characteristics.
As may be observed in two different clusters, the molecules were sorted into groups based on their commonalities ( Figure 6). Molecule 1 is hydroxychloroquine, which presents studies of its activity against the target of SARS-CoV-2. Therefore, the molecules of the cluster in which it is found have the potential to exhibit the same response. Compared to hydroxychloroquine, quinidine, quinine, and tafenoquine were in the second cluster and had the lowest likelihood of acting on the target under study.
According to studies by Da Silva Costa et al. (2018) [15] and Cruz et al. (2018) [16], steric and electrostatic forces influence the structural conformation of molecules, especially those with biological functions, so molecular overlap (overlay) was also carried out for the steric (ste) and electrostatic (ele) fields with the aid of Discovery Studio software (Table 4) using 100% as a parameter.
Hierarchical cluster analysis (HCA) was performed with the help of M software ( State College, PA, USA). A dendrogram of the pharmacophoric chara was generated ( Figure 5), which confirmed the Pearson correlation values, demo a greater approximation of ACC with the Tanimoto index (TI) and a departure HYD and DONN characteristics.  As may be observed in two different clusters, the molecules were sorted into groups based on their commonalities ( Figure 6). Molecule 1 is hydroxychloroquine, which presents studies of its activity against the target of SARS-CoV-2. Therefore, the molecules of the cluster in which it is found have the potential to exhibit the same response. Compared to hydroxychloroquine, quinidine, quinine, and tafenoquine were in the second cluster and had the lowest likelihood of acting on the target under study.  [16], steric and electrostatic forces influence the structural conformation of molecules, especially those with biological functions, so molecular overlap (overlay) was also carried out for the steric (ste) and electrostatic (ele) fields with the aid of Discovery Studio software (Table 4) using 100% as a parameter.   As expected from this investigation, chloroquine provided the best 100ste and 100ele results. For the parameters of 100ste, the three best results of antiparasitics were for quinacrine, quinine, and quinidine, obtaining values of above 80%. Piperaquine, tafenoquine, and amodiaquine performed better for 100ele, with values exceeding 50%. Figure 7 shows the structural comparison of the reference molecule with the other molecules selected for the study. As a result, molecules 9 and 10 differ somewhat in their hydroxyl stereochemistry. In contrast, molecule 11 contains bonds not found in the reference and trifluoride attached to an aromatic ring. These three structures were the ones that significantly differed from the reference due to their chemical characteristics. The main feature common in all structures is the presence of quinoline (C9H7N), a scaffold consistent with the pharmacophores obtained. In addition to having a bonded chlorine atom, molecules 3, 4, and 5 have larger Tanimoto indices and improved structural closeness, as seen in the HCA cluster. Molecule 3 is the closest to the chemical structure of hydroxychloroquine (except for chloroquine).

Virtual Screening
The virtual screening took place on the Pharmit online server using the coordinates obtained from the pharmacophore. Molecular weight (MW), rotatable bonds (RB), partition coefficient (LogS), polar surface area (PSA), number of aromatic rings (AR), hydrogen bond acceptors (HBA), and hydrogen bond donors (HBD) are some of the defined reduction filters that were used to limit the resulting structures to this range of defined values. These filters were selected based on the highest and lowest value presented by all drugs (see Table 5).
The MolPort database was used, and it produced a total of 2000 molecules while accounting for the lowest values of RMSD. All structures showed an RMSD of 0.17, demonstrating the model's effectiveness. Figure 8 illustrates the properties of the pharmacophores for the five compounds obtained from the screening procedure and accepted in this study to illustrate the behavior of the pharmacophore regions. As a result, molecules 9 and 10 differ somewhat in their hydroxyl stereochemistry. In contrast, molecule 11 contains bonds not found in the reference and trifluoride attached to an aromatic ring. These three structures were the ones that significantly differed from the reference due to their chemical characteristics. The main feature common in all structures is the presence of quinoline (C 9 H 7 N), a scaffold consistent with the pharmacophores obtained. In addition to having a bonded chlorine atom, molecules 3, 4, and 5 have larger Tanimoto indices and improved structural closeness, as seen in the HCA cluster. Molecule 3 is the closest to the chemical structure of hydroxychloroquine (except for chloroquine).

Virtual Screening
The virtual screening took place on the Pharmit online server using the coordinates obtained from the pharmacophore. Molecular weight (MW), rotatable bonds (RB), partition coefficient (LogS), polar surface area (PSA), number of aromatic rings (AR), hydrogen bond acceptors (HBA), and hydrogen bond donors (HBD) are some of the defined reduction filters that were used to limit the resulting structures to this range of defined values. These filters were selected based on the highest and lowest value presented by all drugs (see Table 5). The MolPort database was used, and it produced a total of 2000 molecules while accounting for the lowest values of RMSD. All structures showed an RMSD of 0.17, demonstrating the model's effectiveness. Figure 8 illustrates the properties of the pharmacophores for the five compounds obtained from the screening procedure and accepted in this study to illustrate the behavior of the pharmacophore regions.

Prediction of Toxicological and Pharmacokinetic Properties
The 2000 molecules resulting from the screening process were submitted to evaluate their absorption properties, distribution, metabolism, excretion, and toxicity (ADMET). A graph ( Figure 9) was obtained based on the 95% and 99% confidence values for the bloodbrain barrier (BBB), human intestinal absorption (HIA), polar surface area (PSA), and lipid solubility (LogP), which resulted in 1397 molecules for the subsequent phases.

Prediction of Toxicological and Pharmacokinetic Properties
The 2000 molecules resulting from the screening process were submitted to evaluate their absorption properties, distribution, metabolism, excretion, and toxicity (ADMET). A graph (Figure 9) was obtained based on the 95% and 99% confidence values for the blood-brain barrier (BBB), human intestinal absorption (HIA), polar surface area (PSA), and lipid solubility (LogP), which resulted in 1397 molecules for the subsequent phases. Hydroxychloroquine, represented as a star in the graph, was used as a control molecule for pharmacological characteristics. The ellipses show the 95% and 99% confidence-bound space for promising candidates. Those found in this region have features similar to compounds with a capacity of ≥90% to be absorbed, a value of <140 for PSA, and <5 for ALogP98. Only candidates within the 99% confidence limit are selected [17].
The results of the ADMET predictions can be seen in Table 6. The 11 compounds selected together with hydroxychloroquine were evaluated. Ligand_020, ligand_035, ligand_063, and ligand_080, as well as the control drug, showed false negative results for binding to plasma proteins. The others showed false positives and may be highly linked to PPB. This is related to lipophilicity, whereby the higher the lipophilicity, the stronger the link with PPB. In this way, inferences on how these molecules can circulate in the bloodstream can be made. In the same way, PPB properties are related to toxicity that can cause severe unwanted consequences [18]. Table 6. Absorption, distribution, metabolism, excretion, and toxicity (ADMET) predictions for 11 selected compounds and the control (hydroxychloroquine). Hydroxychloroquine, represented as a star in the graph, was used as a control molecule for pharmacological characteristics. The ellipses show the 95% and 99% confidence-bound space for promising candidates. Those found in this region have features similar to compounds with a capacity of ≥90% to be absorbed, a value of <140 for PSA, and <5 for ALogP98. Only candidates within the 99% confidence limit are selected [17].

Molecules
The results of the ADMET predictions can be seen in Table 6. The 11 compounds selected together with hydroxychloroquine were evaluated. Ligand_020, ligand_035, lig-and_063, and ligand_080, as well as the control drug, showed false negative results for binding to plasma proteins. The others showed false positives and may be highly linked to PPB. This is related to lipophilicity, whereby the higher the lipophilicity, the stronger the link with PPB. In this way, inferences on how these molecules can circulate in the bloodstream can be made. In the same way, PPB properties are related to toxicity that can cause severe unwanted consequences [18].
All ligands showed false positive values for hepatotoxicity as with the control drug. This is one of the main drawbacks of synthetic medications, indicating that they can impair liver function or induce liver disorders when used in large doses over an extended period. However, biological studies are needed to prove such data. None of the candidates showed the potential to bind CYP2D6, which is present in 2% of the hepatic CYP content and is responsible for the metabolism of several drugs. Its binding can affect drug metabolism and impair the desired effect [19].
All ligands showed poor or very poor solubility, with emphasis on hydroxychloroquine and ligand_020. Regarding blood-brain barrier penetration, ligand_020 had low penetration, the control drug and ligand_003, ligand_005, and ligand_044 exhibited strong penetration, and the others were just average. Parameters a, b, and c in Table 6 are classified based on a compound with a high binding index (≥90%) with a Bayesian score with a value of −2.209, classified as false positive (true) or false negative (false). The BBB controls the flow of solutes from the blood to the brain and facilitates molecular communication between the central nervous system and other tissues. A high penetration for these candidates is ideal as studies indicate that the SARS-CoV-2 virus can affect brain cells [20].
All molecules in this study showed good intestinal absorption, parameters responsible for drug transport and release. The candidates also underwent in silico toxicity testing using FDA models of male and female rats and mice (see Table 7). Table 6. Absorption, distribution, metabolism, excretion, and toxicity (ADMET) predictions for 11 selected compounds and the control (hydroxychloroquine). None of the candidates showed carcinogenic potential. Hydroxychloroquine had mutagenic potential in the Ames tests, but the study candidates did not, indicating that the chosen compounds may have fewer side effects.

Molecular Docking
Twelve compounds were chosen for the molecular docking investigation after the ADMET evaluation, including hydroxychloroquine, which served as a reference. The simulation was run on the Dockthor server, utilizing coordinates found while determining the hypothetical active site (X = −32.770, Y = 24.940, and Z = 6.690). The molecular target in this study was the structure of the SARS-CoV-2 spike receptor bound with ACE2 (PDB ID: 6M0J). To observe binding affinity values and molecular interactions, see Table 8 vide infra.
Hydroxychloroquine showed a binding affinity of −7.595 kcal·mol −1 to the target and reacted with two residues through hydrogen bonds. The hydrogens of both amines (NH 2 ) of the ARG76 residue interact with the ligand's hydroxyl (-OH). The oxygen atom of ASP73 interacted with the hydroxyl and the amine group of hydroxychloroquine ( Figure 10).
The stability of the ligand within the protein and its binding affinity are significantly influenced by hydrogen bonds, which are of enormous relevance. π-sigma interactions also occur between the aromatic ring with residue THR306 and π-alkyl with the pyridine fragment (C 5 H 5 N) with residues PRO303 and MET365. Amino acids ALA368, ASN304, GLN307, GLY172, GLY336, PHE338, TYR173, and VAL171 interact with the ligand hydrophobically.  The best binding affinity value was observed for ligand_003 with a value of −8.645 kcal·mol −1 . There was only one hydroxyl hydrogen bond for this ligand with residue ASP12 (Figure 11). Other interactions occurred whereby the residue HIS16 formed a π-sulfur interaction with the sulfur present in the structure. The piperidine fragment (C 5 H 11 N) forms a π-alkyl-type interaction with the residues LYS8, PRO371, and VAL75. Additionally, hydrophobic interactions with the residues ARG71, ARG76, ASN15, ASP73, GLN77, GLN78, GLU74, LEU11, LEU123, and LYS85 were observed.
Ligand_033 indicated a binding affinity of −8.303 kcal·mol −1 when bound to the SARS-CoV-2 target structure. This interaction resulted in a hydrogen bond between the amine hydrogen of the ASN15 residue with the pyrimidinone nitrogen (C 4 H 4 N 2 O) present in the ligand. Interaction of the ligand with ARG71 was through a π-cation-type bond ( Figure 12). This was the monopole type (cation), where a large electrostatic force occurs due to the abundance of electrons in the system, being a common interaction in nature. The residue ASP73 formed a carbon-hydrogen interaction with the methyl hydrogens on the ligand. The other interactions were hydrophobic with residues ARG375, ARG76, ASP12, ASP73, GLN77, GLN78, GLU19, GLU74, HIS16, LEU11, LYS8, LYS85, PHE372, PRO371, and TYR173.
For ligand_013, hydrophobic interactions occurred with residues ALA369, ARG76, ASN15, ASP73, GLN370, GLU19, HIS16, PHE372, PRO371, and TYR173. It also formed a hydrogen bond between one of the amine hydrogens of the ARG71 residue with the pyridine nitrogen of the ligand (see Figure 13). A carbon-hydrogen bond occurred between residue ALA369 and the pyrrolidine fragment (C 4 H 9 N) and an π-alkyl type of bond with pyridine. HIS16 and PRO371 residues also made π-alkyl bonds. These interactions occur between the electron cloud and an alkyl group. HIS16 also formed a stacked π-π interaction. These stacked π-π interactions are of the dipole-dipole type occurring due to the electrostatic attraction of the electrons of the aromatic ring. They play an important role in the stability of proteins. GLU19 makes a π-anion-type bond, occurring mainly in electron-deficient aromatic rings. This interaction is rare in biological systems as amino acids tend to repel nearby negative charges. The stability of the ligand within the protein and its binding affinity are significantly influenced by hydrogen bonds, which are of enormous relevance. π-sigma interactions also occur between the aromatic ring with residue THR306 and π-alkyl with the pyridine fragment (C5H5N) with residues PRO303 and MET365. Amino acids ALA368, ASN304, GLN307, GLY172, GLY336, PHE338, TYR173, and VAL171 interact with the ligand hydrophobically.
The best binding affinity value was observed for ligand_003 with a value of −8.645 kcal.mol −1 . There was only one hydroxyl hydrogen bond for this ligand with residue ASP12 (Figure 11). Other interactions occurred whereby the residue HIS16 formed a πsulfur interaction with the sulfur present in the structure. The piperidine fragment (C5H11N) forms a π-alkyl-type interaction with the residues LYS8, PRO371, and VAL75. Additionally, hydrophobic interactions with the residues ARG71, ARG76, ASN15, ASP73, GLN77, GLN78, GLU74, LEU11, LEU123, and LYS85 were observed.
Ligand_033 indicated a binding affinity of −8.303 kcal.mol −1 when bound to the SARS-CoV-2 target structure. This interaction resulted in a hydrogen bond between the amine hydrogen of the ASN15 residue with the pyrimidinone nitrogen (C4H4N2O) present in the ligand. Interaction of the ligand with ARG71 was through a π-cation-type bond ( Figure  12). This was the monopole type (cation), where a large electrostatic force occurs due to the abundance of electrons in the system, being a common interaction in nature. The residue ASP73 formed a carbon-hydrogen interaction with the methyl hydrogens on the ligand. The other interactions were hydrophobic with residues ARG375, ARG76, ASP12, ASP73, GLN77, GLN78, GLU19, GLU74, HIS16, LEU11, LYS8, LYS85, PHE372, PRO371, and TYR173. Ligand_044 showed a binding affinity of −7.749 kcal·mol −1 and formed two hydrogen bonds between the oxygen of the ligand and hydrogen of the amine present in residue ARG71 and residue LYS85 ( Figure 14). Hydrophobic interactions occurred with 16 amino acid residues, ALA368, ALA369, ARG375, ARG76, ASN15, ASP73, GLN370, GLN77, GLU19, GLU74, GLY84, HIS16, LYS85, PHE372, PRO371, and TYR173. The ALA369 and LYS85 residues also form π-alkyl interactions with the azepane fragment (C 6 H 13 N) and the pyrimidines. The HIS16 residue forms a π-π-type interaction stacked with the aromatic ring, while the amino acid TYR173 forms a bond of the carbon-hydrogen type. There are also halogenated interactions of fluorine with the amino acids ASN15 and GLU19. The homeostasis of various physiological processes is directly influenced by fluorine. Halogens in the ligands may improve drug binding affinity and selectivity by acting as an electron acceptor and directly participating in the interaction with biological receptors.
The binding affinity of ligand_080 was −7.690 kcal·mol −1 . The ligand generated 14 hydrophobic interactions with the Spike-ACE2 target with residues ARG71, ARG76, ASN15, ASP73, GLN78, GLU74, GLY84, HIS16, LEU11, LYS8, LYS85, PRO371, THR83, and VAL75. Three hydrogen bonds were formed between the oxygen atom of ASP12 and the nitrogen in the ligand structure, the amine nitrogen of GLN77 and ligand oxygen, and one with residue GLN78 (Figure 15).     bond with pyridine. HIS16 and PRO371 residues also made π-alkyl bonds. These interactions occur between the electron cloud and an alkyl group. HIS16 also formed a stacked π-π interaction. These stacked π-π interactions are of the dipole-dipole type occurring due to the electrostatic attraction of the electrons of the aromatic ring. They play an important role in the stability of proteins. GLU19 makes a π-anion-type bond, occurring mainly in electron-deficient aromatic rings. This interaction is rare in biological systems as amino acids tend to repel nearby negative charges. Ligand_044 showed a binding affinity of −7.749 kcal.mol −1 and formed two hydrogen bonds between the oxygen of the ligand and hydrogen of the amine present in residue ARG71 and residue LYS85 ( Figure 14). Hydrophobic interactions occurred with 16 amino acid residues, ALA368, ALA369, ARG375, ARG76, ASN15, ASP73, GLN370, GLN77, GLU19, GLU74, GLY84, HIS16, LYS85, PHE372, PRO371, and TYR173. The ALA369 and LYS85 residues also form π-alkyl interactions with the azepane fragment (C6H13N) and the pyrimidines. The HIS16 residue forms a π-π-type interaction stacked with the aromatic ring, while the amino acid TYR173 forms a bond of the carbon-hydrogen type. There are also halogenated interactions of fluorine with the amino acids ASN15 and GLU19. The homeostasis of various physiological processes is directly influenced by fluorine. Halogens in the ligands may improve drug binding affinity and selectivity by acting as an electron acceptor and directly participating in the interaction with biological receptors.  By examining the docking results for all 12 ligands and hydroxychloroquine, it was feasible to confirm that more hydrophobic interactions (34 residues) than hydrogen bonds (11 residues) were used to make contact between the ligands and the Spike-ACE2 protein (see Figure 16).  By examining the docking results for all 12 ligands and hydroxychloroquine, it was feasible to confirm that more hydrophobic interactions (34 residues) than hydrogen bonds (11 residues) were used to make contact between the ligands and the Spike-ACE2 protein (see Figure 16). By examining the docking results for all 12 ligands and hydroxychloroquine, it was feasible to confirm that more hydrophobic interactions (34 residues) than hydrogen bonds (11 residues) were used to make contact between the ligands and the Spike-ACE2 protein (see Figure 16).  Interacting residues are predominant between the region of α-helix A and β-sheet A of ACE2, and the region of β-sheet B of the spike protein ( Figure 17). The most frequent residue was PRO371, which was present in β-sheet A in 77% of the cases. Next in frequency were residues HIS16 from the α-helix and TYR173 found in β-sheet A, both of which were detected in 69% of the instances. Interacting residues are predominant between the region of α-helix A and β-sheet A of ACE2, and the region of β-sheet B of the spike protein ( Figure 17). The most frequent residue was PRO371, which was present in β-sheet A in 77% of the cases. Next in frequency were residues HIS16 from the α-helix and TYR173 found in β-sheet A, both of which were detected in 69% of the instances. When all interactions of each ligand are taken into account, the residue with the highest frequency for hydrogen bonds was ARG71 with 31% in β-sheet B, which is favorable to what is anticipated because this region comes from the spike bond, followed by residues ALA368 (β-sheet A), ARG375 (β-sheet A), ASN15 (α-helix A), and ASP12 (αhelix A) with 15%. All amino acid residues that established hydrogen bonds also interacted through hydrophobic interactions.
Arginine, a positively charged guanidinium group with pH 7.0, is a residue responsible for catalyzing several enzymes. Its main action is found in the synthesis of nitric oxide. The interaction of the ligand with this residue may favor the inhibition of the spike protein, as it is responsible for decreasing the bioavailability of nitric oxide, and responsible for vasodilation and blood pressure control.

Synthetic Accessibility Prediction and Similarity Analysis
Following the ADMET examination, the 99 compounds were put through the synthetic accessibility (SA) test, which filtered out the molecules that were not likely to be synthesized. The AMBIT-AS program was used to predict SA, and its algorithm computes scores ranging from 0 to 100 based on several criteria relating to the chemical properties of the structures, where a molecule with great ease of synthesis has a value of 100 [21]. Hydroxychloroquine was used as a standard and had an SA value of 69.792. Here, we only selected the molecules that achieved values greater than this reference molecule, indicating that they were highly likely to be synthesized (see Table 9).  When all interactions of each ligand are taken into account, the residue with the highest frequency for hydrogen bonds was ARG71 with 31% in β-sheet B, which is favorable to what is anticipated because this region comes from the spike bond, followed by residues ALA368 (β-sheet A), ARG375 (β-sheet A), ASN15 (α-helix A), and ASP12 (α-helix A) with 15%. All amino acid residues that established hydrogen bonds also interacted through hydrophobic interactions.
Arginine, a positively charged guanidinium group with pH 7.0, is a residue responsible for catalyzing several enzymes. Its main action is found in the synthesis of nitric oxide. The interaction of the ligand with this residue may favor the inhibition of the spike protein, as it is responsible for decreasing the bioavailability of nitric oxide, and responsible for vasodilation and blood pressure control.

Synthetic Accessibility Prediction and Similarity Analysis
Following the ADMET examination, the 99 compounds were put through the synthetic accessibility (SA) test, which filtered out the molecules that were not likely to be synthesized. The AMBIT-AS program was used to predict SA, and its algorithm computes scores ranging from 0 to 100 based on several criteria relating to the chemical properties of the structures, where a molecule with great ease of synthesis has a value of 100 [21]. Hydroxychloroquine was used as a standard and had an SA value of 69.792. Here, we only selected the molecules that achieved values greater than this reference molecule, indicating that they were highly likely to be synthesized (see Table 9).
After making accessibility predictions, the Chemmine Tools server was used to conduct similarity research. The structures were divided into clusters according to their individual physical-chemical and structural features. They were also divided into clusters, and five molecular structures-one for each cluster-and hydroxychloroquine was chosen to proceed with this study. The selection was based on the best binding affinity values, taking into account that there is a chance that they will have similar biological activity in in vitro experiments if they are in the same cluster (see Figure 18). Ultimately, 12 interesting compounds were chosen for further study, and molecular docking studies were run on each.

Lipophilicity and Water Solubility
The SwissADME server was employed to ascertain the selected molecules' lipophilicity and water solubility characteristics. These properties have an impact on the kinetics of how drug candidates act, making them crucial in the process of identifying novel medications. The octanol-water partition coefficient (LogPo/w), which results from the concentration of a molecule in its neutral form in the organic and aqueous phase, is used as a benchmark for lipophilicity [18]. Lipophilicity affects the solubility and permeability properties of membranes in a way that can change the molecules' ADMET profiles and directly reflect the molecular forces active between the two phases. For novel compounds, LogP < 5 is a desirable number [18]. After making accessibility predictions, the Chemmine Tools server was used to conduct similarity research. The structures were divided into clusters according to their individual physical-chemical and structural features. They were also divided into clusters, and five molecular structures-one for each cluster-and hydroxychloroquine was chosen to proceed with this study. The selection was based on the best binding affinity values, taking into account that there is a chance that they will have similar biological activity in in vitro experiments if they are in the same cluster (see Figure 18). Ultimately, 12 interesting compounds were chosen for further study, and molecular docking studies were run on each.

Lipophilicity and Water Solubility
The SwissADME server was employed to ascertain the selected molecules' lipophilicity and water solubility characteristics. These properties have an impact on the kinetics of how drug candidates act, making them crucial in the process of identifying novel medications. The octanol-water partition coefficient (LogPo/w), which results from the concentration of a molecule in its neutral form in the organic and aqueous phase, is The SwissADME server allows the investigation of five distinct values, XLOGP3 which uses a database of molecules as a basis for identifying LogP values, using 87 different types of atoms and two groups of correction factors, WLOGP that uses a fragmented method proposed by Wildman and Crippe through atomic contributions, MLOGP that uses 13 descriptors of hydrophobic, hydrophilic, unsaturated bond, and other atoms based on the linear relationship to determine their values, SILICOS-IT which uses seven descriptors in its approach in a hybrid system based on the FILTER-IT software source code, and iLOGP that uses the energies of free water and n-octanol solvations using a generalized Born equation and solvent accessible surface area (GB/AS) approach.
The five selected structures and hydroxychloroquine were passed through the Swis-sADME server to obtain their lipophilicity values and an average of the results (Table 10 and Figure 19). As suggested by Lipinski for compounds with optimal lipophilicity, all molecules displayed values below 5, making them desirable. Ligand_044 presented a higher average than the control drug, which could be associated with the presence of the halogen (F) in its structure. Hydroxychloroquine, the standard, also has a halogen (Cl) registered as the second highest value. The third highest result was ligand_003, possibly due to the double bond with sulfur (S). The values obtained in these studies ranged from 2.58 to 4.22, characterizing them as highly lipophilic, meeting the criteria required for drug candidates.
Water solubility (LogS) is a crucial investigational factor since it relates to the route of delivery of drug candidates, particularly when addressing oral and parenteral routes, due to the concentration required for the candidate's circulation in the distribution system. The usage of molecules with low solubility has proven problematic since their insoluble classification affects how they behave in the circulatory system. Many non-polar alkene and alkyne bonds in the structures affect the compounds' solubility since they prevent the compounds from dissolving in water. The SwissADME server allows the investigation of five distinct values, XLOGP3 which uses a database of molecules as a basis for identifying LogP values, using 87 different types of atoms and two groups of correction factors, WLOGP that uses a fragmented method proposed by Wildman and Crippe through atomic contributions MLOGP that uses 13 descriptors of hydrophobic, hydrophilic, unsaturated bond, and other atoms based on the linear relationship to determine their values, SILICOS-IT which uses seven descriptors in its approach in a hybrid system based on the FILTER-IT software source code, and iLOGP that uses the energies of free water and n-octanol solvations using a generalized Born equation and solvent accessible surface area (GB/AS) approach.
The five selected structures and hydroxychloroquine were passed through the SwissADME server to obtain their lipophilicity values and an average of the results (Table  10 and Figure 19). As suggested by Lipinski for compounds with optimal lipophilicity, al molecules displayed values below 5, making them desirable. Ligand_044 presented a higher average than the control drug, which could be associated with the presence of the halogen (F) in its structure. Hydroxychloroquine, the standard, also has a halogen (Cl registered as the second highest value. The third highest result was ligand_003, possibly due to the double bond with sulfur (S). The values obtained in these studies ranged from 2.58 to 4.22, characterizing them as highly lipophilic, meeting the criteria required for drug candidates.  Water solubility (LogS) is a crucial investigational factor since it relates to the route of delivery of drug candidates, particularly when addressing oral and parenteral routes due to the concentration required for the candidate's circulation in the distribution The molecules are more soluble in water when a hydroxyl group is present, and a ketone functional group has a similar effect because of the potential for molecular interactions. With the help of SwissADME, three methods-the ESOL method, the Ali method, and the SILICOS-IT technique-were designed for this study to determine water solubility [22].
In terms of solubility, ligand_044 and hydroxychloroquine produced the best results. According to the SILICOS-IT parameter, ligand_044 (−7.46) and hydroxychloroquine (−6.35) have low solubility. A moderate solubility displays a LogS between −4 and −6, while a high solubility displays a LogS from −2 to −4. Except for ligand_013 and ligand_044, all others showed moderate water solubility. However, they are still potential candidates for oral or parenteral administration (Table 11 and Figure 20). The molecules are more soluble in water when a hydroxyl group is present, and a ketone functional group has a similar effect because of the potential for molecular interactions. With the help of SwissADME, three methods-the ESOL method, the Ali method, and the SILICOS-IT technique-were designed for this study to determine water solubility [22].
In terms of solubility, ligand_044 and hydroxychloroquine produced the best results. According to the SILICOS-IT parameter, ligand_044 (−7.46) and hydroxychloroquine (-6.35) have low solubility. A moderate solubility displays a LogS between −4 and −6, while a high solubility displays a LogS from −2 to −4. Except for ligand_013 and ligand_044, all others showed moderate water solubility. However, they are still potential candidates for oral or parenteral administration (Table 11 and Figure 20).

Molecular Dynamics Studies
To understand the stability and binding affinity of ligands complexed in the spike protein binding site, we performed an MD simulation. Protein conformational changes disclosed its function, and MD is crucial to comprehending biological activity. After MD simulation, trajectory data were analyzed in terms of RMSD, RMSF, radius of gyration, VMD, and R-code. RMSD based on Cα of spike protein with its complexes is plotted in Figure 21.
The RMSD plot shows very stable curves that indicate the system's stability during the trajectory. The average RMSD of all ligands has similar curves compared with the control compound (hydroxychloroquine) and ranges from 2-4 Å. The effect of conformational change is more (2-4 Å) in the control compound, ligand_003-Spike-ACE2 complex, and ligand_080-Spike-ACE2 complex. Ligand_033-Spike-ACE2 complex was moderate (2-3.5 Å) and ligand_013-Spike-ACE2 complex and ligand_044-Spike-ACE2 complex (2-3 Å) were lowest. Low variation indicates that ligands are still interacting within the binding cavity of the spike protein during the simulation period, despite minor conformational changes.

Molecular Dynamics Studies
To understand the stability and binding affinity of ligands complexed in the spike protein binding site, we performed an MD simulation. Protein conformational changes disclosed its function, and MD is crucial to comprehending biological activity. After MD simulation, trajectory data were analyzed in terms of RMSD, RMSF, radius of gyration, VMD, and R-code. RMSD based on Cα of spike protein with its complexes is plotted in Figure 21.
shows more flexibility due to the high degree of freedom. The average fluctuation 1-2 Å at the central residues, while the loop region shows fluctuation up to 5 Å that the system's overall stability. Ligand_033-Spike-ACE2 complex shows a little fluct at residual position 400 during simulation, but that residue does not share a binding region, so its variation does not affect binding strength of the complex.  The RMSD plot shows very stable curves that indicate the system's stability during the trajectory. The average RMSD of all ligands has similar curves compared with the control compound (hydroxychloroquine) and ranges from 2-4 Å. The effect of conformational change is more (2-4 Å) in the control compound, ligand_003-Spike-ACE2 complex, and ligand_080-Spike-ACE2 complex. Ligand_033-Spike-ACE2 complex was moderate (2-3.5 Å) and ligand_013-Spike-ACE2 complex and ligand_044-Spike-ACE2 complex (2-3 Å) were lowest. Low variation indicates that ligands are still interacting within the binding cavity of the spike protein during the simulation period, despite minor conformational changes.
The Cα-based RMSF was observed, as shown in Figure 22. The overall RMSF curve of all the complexes is stable from the central region, especially the binding site region, while fluctuating at the loop region, including residues from 600 to 800. The loop region shows more flexibility due to the high degree of freedom. The average fluctuation is from 1-2 Å at the central residues, while the loop region shows fluctuation up to 5 Å that shows the system's overall stability. Ligand_033-Spike-ACE2 complex shows a little fluctuation at residual position 400 during simulation, but that residue does not share a binding cavity region, so its variation does not affect binding strength of the complex. The compactness of the protein is shown by the radius of gyration, as indicated in Figure 23. Based on mean values that range from 49.3 to 49.6 Å, it can be concluded that the protein is still compact and that no significant modifications are seen when spike protein is in the presence of its ligands. The constant mean value shows that protein is stable and there are no major conformational changes in spike protein. According to RMSD, RMSF, and gyration radius, the ligands could stabilize the spike protein binding site during the MD trajectories. The compactness of the protein is shown by the radius of gyration, as indicated in Figure 23. Based on mean values that range from 49.3 to 49.6 Å, it can be concluded that the protein is still compact and that no significant modifications are seen when spike protein is in the presence of its ligands. The constant mean value shows that protein is stable and there are no major conformational changes in spike protein. According to RMSD, RMSF, and gyration radius, the ligands could stabilize the spike protein binding site during the MD trajectories.

Binding Free Energy Calculation
The binding energy of ligand_003, 013, 033, 044, and 080 in the spike protein binding site was calculated by MMPBSA methods, and obtained results are summarized in Table 12. The ∆G total indicated that all the complexes had excellent simulation-based binding energy values that were near to those of the control compound. The MMPBSA result ranks the compounds ligand_44 < ligand_33 < hydroxychloroquine (control) < ligand_03, ligand_80 < ligand_013.
Moreover, different components of the energies, vdW, electrostatics, polar solvation energy (EPB), and non-polar solvation energy (ENPOLAR), all show good binding strength. In control, compound electrostatics is a major energy contribution, whereas in other compounds, vdW is dominant over other energy terms.
The fact that the polar solvation energy (EPB) remained notably unfavorable also indicates that vdW plays a key role in the binding of ligands to spike protein. Ligand_044 shows a stronger and more significant free energy value than other complexes and has the highest vdW, electrostatic, and EPB values. While all compounds had positive binding interactions and significant binding energy values in the protein binding site, ligand_013 had the lowest ∆G total due to the lowest electrostatic contribution.

Binding Free Energy Calculation
The binding energy of ligand_003, 013, 033, 044, and 080 in the spike protein binding site was calculated by MMPBSA methods, and obtained results are summarized in Table  12. The ∆Gtotal indicated that all the complexes had excellent simulation-based binding energy values that were near to those of the control compound. The MMPBSA result ranks the compounds ligand_44 < ligand_33 < hydroxychloroquine (control) < ligand_03, ligand_80 < ligand_013.

Determination of the Theoretical Mean Inhibitory Concentration (IC 50 )
The mean inhibitory concentration (IC 50 ) is a measure used to assess the ability of a compound to inhibit half of the biological activity in a specific target from an initial amount [23,24]. The IC 50 measurement value has always been one of the limitations of theoretical and in silico approaches. In this work, we used a mathematical equation proposed by Hopkins et al. (2014) [25], in which the initial formula was deduced and rearranged to estimate the pIC 50 value (Table 13). As is already known for pIC 50 values, the higher the value obtained, the more active the compound studied, while the IC 50 value is the opposite, i.e., the lower its value, the better the activity. With a predicted IC 50 value of 0.459 µM, ligand_080 has the potential to be the most active. The values discovered were 0.586 µM and 0.663 µM for ligand_013 and ligand_033, respectively. For ligand_003, the highest value was 2.371 µM. The value of hydroxychloroquine was inferior to all the predicted values for the ligands.

Quantum Chemical Calculations
Chemical descriptors play a crucial role in the chemistry and pharmacology of the interaction between the ligand and the macromolecule. The HOMO and LUMO energies and their gap show the charge transfer potential that affects the molecule's bioactivation and the electrostatic potential connected to the binding of ligands to the surface of the macromolecule [26].
The five candidate compounds were subjected to chemical molecular optimization calculations using the Gaussian 09 program. As a result, the highest occupied molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs) of each structure, as well as their gap values, were determined for each compound ( Figure 24). The HOMOs and LUMOs are of great importance in the study of the reactivity and stability of molecules, where the HOMO energy demonstrates the ability of a molecule to donate its electrons and the LUMO energy to accept electrons.
The HOMOs are predominant mainly over the pyridine and pyrimidine fragments of the structures, except ligand_080 where this fragment did not present a molecular orbital, a functional group known to perform hydrophobic interactions in its aromatic rings, which was observed in molecular docking studies. The presence of these orbitals explains the π-anion bond present in the aromatic ring of ligand_013. For ligands that contain sulfur, it is possible to see the orbital density on this atom favoring the π-sulfur bond as observed in docking studies, as well as the fluorine atoms that formed the halogen bonds with the spike target. The oxygen atoms also presented HOMOs except for the hydroxyl in ligand_003 and the ketone in ligand_080. The eigenvalues were all very similar, showing that their electron-donating potential is similar.
The HOMOs and LUMOs are localized on the pyridine and pyrimidine fragments with the exception of some fragments that mainly have rings linked to fluorine. For the LUMO values, it is possible to see a difference where ligand_033 had the lowest value of −1.638 eV and ligand_080 the highest of −2.699 eV. The orbitals on ligand_033 favored the π-cation bond, which is explained by the density of the orbital on the aromatic ring next to the oxygen of the ether. calculations using the Gaussian 09 program. As a result, the highest occupied molecular orbitals (HOMOs) and lowest unoccupied molecular orbitals (LUMOs) of each structure, as well as their gap values, were determined for each compound (Figure 24). The HOMOs and LUMOs are of great importance in the study of the reactivity and stability of molecules, where the HOMO energy demonstrates the ability of a molecule to donate its electrons and the LUMO energy to accept electrons. The HOMOs are predominant mainly over the pyridine and pyrimidine fragments of the structures, except ligand_080 where this fragment did not present a molecular orbital, a functional group known to perform hydrophobic interactions in its aromatic rings, which was observed in molecular docking studies. The presence of these orbitals explains the π-anion bond present in the aromatic ring of ligand_013. For ligands that contain sulfur, it is possible to see the orbital density on this atom favoring the π-sulfur bond as observed in docking studies, as well as the fluorine atoms that formed the halogen bonds with the spike target. The oxygen atoms also presented HOMOs except for the hydroxyl in ligand_003 and the ketone in ligand_080. The eigenvalues were all very similar, showing that their electron-donating potential is similar.
The HOMOs and LUMOs are localized on the pyridine and pyrimidine fragments with the exception of some fragments that mainly have rings linked to fluorine. For the LUMO values, it is possible to see a difference where ligand_033 had the lowest value of −1.638 eV and ligand_080 the highest of −2.699 eV. The orbitals on ligand_033 favored the π-cation bond, which is explained by the density of the orbital on the aromatic ring next to the oxygen of the ether.
These orbitals are of great importance in biological interactions, mainly influencing the binding affinity of ligands with targets through interactions, as seen in molecular docking. The HOMO and LUMO data were also important in the calculation of the descriptors gap energy (∆E), ionization potential (IP), electron affinity (EA), These orbitals are of great importance in biological interactions, mainly influencing the binding affinity of ligands with targets through interactions, as seen in molecular docking. The HOMO and LUMO data were also important in the calculation of the descriptors gap energy (∆E), ionization potential (IP), electron affinity (EA), electronegativity (χ), chemical potential (µ), chemical hardness (η), softness (σ), and electrophilicity (ω), see Table 14.
The values of ∆E are related to the chemical reactivity of the compound and are obtained from the difference in the energy of HOMO and LUMO, where a smaller gap indicates the low stability of the structure. According to the data collected, ligand_013 has the best stability among the potential ligands, attaining a value of 4.380 ev. The others presented values that vary between 3.449 eV and 4.434 eV and this demonstrates that all candidates have excellent chemical stability, favoring the stability of the ligands within the spike target substrate.
The IP is related to the ability of the structures studied to donate electrons. Since it relates to the transfer of electrons present in the structure, it is a parameter that may be used to determine the antioxidant potential of the candidates evaluated. All compounds had close values, with the highest value of 6.334 for ligand_013 and the lowest value of 6.071 for ligand_033. This demonstrates that all compounds have an antioxidant power and are great candidates for anti-SARS-CoV-2 action. EA is related to the ability of the ligands to accept electrons. Its relationship to antioxidant activity is unknown. However, it has a characteristic that is the opposite of IP.
Its connection to the potential for reduction and its capacity to neutralize free radicals is apparent. All candidates had a low value of EA, as expected. The lowest value observed was for ligand_033 of 1.638 eV and the highest for ligand_080 (2.699 eV), demonstrating the low potential for accepting electrons and confirming its antioxidant power.
The µ is related to the energetic change of the system through the electronic trend influencing the molecule's electron density with the target substrate. The highest observed value was −4.423 eV for ligand_080, and the lowest was −3.855 eV for ligand_033. The values were approximate for all ligands showing a nucleophilic (electron donor) and electrophilic (electron acceptor) character indicating their charge transfer capacity.
The η is related to the molecular stability of compounds following the maximum hardness principle (MHP), which indicates that hard molecules will be less reactive than soft molecules due to resistance in charge transfer. All compounds showed similar values indicating that the molecules have low hardness and are reactive. The acceptance of loads and measures of the degree of molecular activity characterize the parameter σ. The studied ligands showed a high smoothness.
The value of ω relates to the direction in which the charge transfer occurs and its electrophilic potential, that is, the molecule's stability when acquiring an extra electron density. The molecules studied can thus be classified as strong electrophiles.
The values of χ are related to the attraction of electrons in the molecule. Another piece of information that can be observed is the molecule's electrostatic potential, which can be observed in Figure 25. The ligands were shown to be electronegative, which provides the exchange of electrons with the substrate, favoring electrostatic attraction and increasing the force of interaction.  The electrophilic and nucleophilic areas of each ligand's structure can be observed using the molecular electrostatic potential, with the most electronegative region being represented by red and the least electronegative part by blue. These areas are connected to the area where the majority of a compound's molecular interactions occur. In ligand_003, the most electronegative region is close to oxygen and sulfur, the same region that presented the hydrogen bond with the spike protein. The least electronegative region was at the ends of the molecule in the aromatic ring and pyridine. For ligand_013, there was also a higher electronegative concentration of oxygen; as observed in the molecular docking in this region, hydrogen bonds to the spike protein occurred.
Along with the other ligands, ligand_033 displayed the electronegative region on the The electrophilic and nucleophilic areas of each ligand's structure can be observed using the molecular electrostatic potential, with the most electronegative region being represented by red and the least electronegative part by blue. These areas are connected to the area where the majority of a compound's molecular interactions occur. In ligand_003, the most electronegative region is close to oxygen and sulfur, the same region that presented the hydrogen bond with the spike protein. The least electronegative region was at the ends of the molecule in the aromatic ring and pyridine. For ligand_013, there was also a higher electronegative concentration of oxygen; as observed in the molecular docking in this region, hydrogen bonds to the spike protein occurred.
Along with the other ligands, ligand_033 displayed the electronegative region on the structure's two oxygens and the opposite region on the pyridine and pyrimidine nitrogens, the region in which the hydrogen bonds took place. While ligand_044 displayed an electronegative density in the nitrogen and oxygen regions, both sites for hydrogen bonding, fluorine did not exhibit a high electronegativity for this ligand. For ligand_080, the most electronegative region is found on the oxygens and fluorine, as expected, the regions where hydrogen bonds occurred in molecular docking.

Determination of Theoretical Active Site
With the use of the FTMap server, identifying the potential region of biological activity was investigated (https://ftmap.bu.edu/, accessed on 1 January 2022) [27], whereby regions with significant contributions to the energy of the ligand interaction with hot spots in the macromolecule were identified. Small organic compounds are used to map the target protein by acting as probes on its surface and locating critical regions of interest.
Hot spots have already been studied as one of the ways to identify active sites, especially in structures that do not have co-crystalized ligands. These regions act as indicators of where a drug is likely to bind. The protein structure's cavities were further examined using Molegro software [28][29][30], bearing in mind that studies suggest that the cavities are situated near active sites.

Selection of Antiparasitics and Similarity of Tanimoto
Nine antiparasitic drugs were selected from the DrugBank server (https://go.drugbank. com/, accessed on 1 March 2022) [31], evaluated for their structural similarity with hydroxychloroquine, and their inhibition of SARS-CoV-2 spike protein was determined. The BindingDB server platform (http://bindingdb.org/bind/index.jsp, accessed on 5 March 2022) [32] was used to calculate the Tanimoto index as per the equation vide infra. The software utilizes merged similarity scoring, where the numbers of bits in x and y are set to 1.
The outcome rates a molecule's similarity to another that exhibits similar properties. Therefore, the likelihood that this molecule will be active increases as M(x) increases.

Hierarchical Cluster Analysis (HCA) and Molecular Overlay
The pharmacophoric descriptors were employed in the adopted methodology [35] to assess their link to the Tanimoto index values, observing the significance of each factor for similarity through the Euclidean distance correlation with Minitab 19 software. The Euclidean distance, one of the most used models, is used to carry out clustering through models to identify commonalities between two points. Since d ab is the coordinate of the point in the r dimension, the Euclidean distance between points a and b can be calculated by the formula: 1 2 where p denotes how many dimensions there are in space. The distance between points a and b in the coordinate system can therefore be calculated by applying the equation mentioned above. Here, the hierarchical cluster analysis (HCA) was used to verify the degree of similarity between the molecules studied.

Virtual Screening and Selection of Inhibitory Compounds
The 2000 molecules closest to the adopted models were chosen using the pharmacophoric characteristics generated, which included molecular weight (MW), rotatable bonds (RB), water solubility (LogS), polar surface area (PSA), number of aromatic rings (AR), hydrogen bond acceptors (HBA), and hydrogen bond donors (HBD). This virtual screening was carried out in the MolPort database using the Pharmit server (http: //pharmit.csb.pitt.edu/, accessed on 7 May 2022) [36].

Prediction of Toxicological and Pharmacokinetic Properties
The structures obtained from the virtual screening were subjected to in silico predictions of pharmacokinetic and toxicological features, or absorption, distribution, metabolism, excretion, and toxicity (ADMET). The computations utilize physicochemical factors, drug similarity, and pharmacokinetic profiles to calculate prediction values [37][38][39]. The drugs hydroxychloroquine and chloroquine were used as a benchmark.

Molecular Docking
This study utilized the protein structure of the SARS-CoV-2 spike receptor binding domain co-crystalized with ACE2 (Spike-ACE2) (PDB ID: 6M0J) [40] obtained from the Protein Data Bank database (https://www.rcsb.org/, accessed on 25 May 2022). The UCSF Chimera software [41] was used to remove water molecules and other residues resulting from crystallography that could interfere with the ligand-macromolecule interaction. The molecules with the best pharmacokinetic and toxicological parameters were selected for molecular docking, and hydroxychloroquine and chloroquine were used as controls.
The protein was prepared in APBS-PDB2PQR (https://server.poissonboltzmann.org/, accessed on 25 May 2022) [42], maintaining neutral pH and hence simulating the pH of the organism. The PARSE force field was used to correct the amino acid chains and adjust the conformation. Molecular docking simulations were performed through the DockThor server (https://www.dockthor.lncc.br/v2/, accessed on 26 May 2022) [43]. The grid was selected based on the x, y, and z coordinates of the active site obtained by determining the hot spots. A 20 × 20 × 20 cm cubic box was used, and the other standard parameters (number of evaluations 500,000, population size 750, and 12 runs) were used to analyze the conformations, interaction of molecules with protein amino acids, and binding energy.

Synthetic Accessibility Prediction and Similarity Analysis
Synthetic accessibility (SA) prediction was performed for the molecules predicted to have the best pharmacokinetic and toxicological properties through the AMBIT-AS software (http://ambit.sourceforge.net/reactor.html, accessed on 7 June 2022). This server analyzes the structural and topological characteristics of the molecules studied based on stereochemistry. The algorithm works with a score ranging from 0 to 100, where 100 symbolizes a molecule that is easily synthesized.
The similarity analysis of the molecules was evaluated with the Chemmine Tools server (https://chemminetools.ucr.edu/, accessed on 7 June 2022) through hierarchical clustering [44]. The Tanimoto index was used to calculate the similarity through the atomic descriptors through a matrix of the unique characteristics of each structure [45,46].

Molecular Dynamics
Molecular dynamic (MD) simulation was performed for 100ns using NAMD software [48,49] in order to study the stability and binding energy of the spike protein model complexed with ligand_4, ligand_33, ligand_03, ligand_80, ligand_013. The input files were prepared using antechamber and t-leap modules of Amber 14 tools [50] and then, a 50ns MD simulation was run for five complexes along with the control compound. The LEAP program of Amber 14 [51] was employed to produce the force field, coordinate, and topology information of the complexes. Generalized Amber Force Field (GAFF) [52] was used to generate the ligand parameters, while ff14SB force field [53] was employed for protein parameters. The system was solvated by the TIP3P water model and neutralized by adding Na+ ions. The final system comprises water, ligand, and protein complexes.
To maintain the protein's constraint at its mean position, the systems were minimized using the steepest descent minimization methodology. After that, the system's temperature gradually increased by 300K under ensemble conditions (NPT). In the equilibrium phase, electrostatic interactions were enumerated using the particle-mesh Ewald algorithm [54] with 12 Å cutoff distance and periodic boundary conditions employed to calculate the forces of atoms, whereas all bonded interactions, such as hydrogen bonds, were constrained by the SHAKE algorithm [55]. Langevin [56] coupling was employed at constant temperature. The trajectories data were computed, and a snapshot taken at each 20ps time step.
MD trajectory was analyzed using molecular dynamic software (VMD), R-program, and Pymol. The stability and binding energy of the complexes were calculated by RMSD, RMSF, and radius of gyration. RMSD data provide average motion between coordinates.

Binding Free Energy Calculation
The free energy calculations were performed using the molecular mechanics energies combined with Poisson-Boltzmann (MM-PBSA) [57][58][59]. The free energy was calculated as follows: where ∆G bind is the free energy of the complex, resulting from the sum of the molecular mechanics' energy (∆E MM ), desolvation free energy (∆G solv ), and entropy (−T∆S).
The energy of molecular gas phase mechanics (∆E MM ) can be described by the sum of the internal energy contributions (∆E internal ); the sum of the connection, angle, and dihedral energies; electrostatic contributions (∆E electrostatic ); and van der Waals terms (∆E vdW ).
The desolvation free energy (∆G solv ) is the sum of the polar (∆G GB ) and non-polar (∆G nonpol ) contributions. The polar desolvation term was calculated using the implicit generalized Born (GB) approach. g_mmpbsa was employed to investigate the binding free energy of the selected complex at the binding site of the spike protein. The binding free energy was decomposed into relative free energy of solvated complex (protein and ligand) and discrete receptor and ligand components given by the equation: For each free energy, it is a summarization of different molecular mechanics energy including polar and non-polar solvation energy, electrostatic, and van der Waal's contributions.

Determination of Theoretical Mean Inhibitory Concentration (IC 50 )
After the binding affinity values were determined, they were utilized to calculate the ligand efficiency (LE) value [23,24], where hydroxychloroquine was employed as the study's point of reference, acting as a filter in the assessment of energies. This procedure was applied to the selected molecules; the efficiency was evaluated using the equation: where ∆G is the energy obtained from molecular docking and N is the number of nonhydrogen atoms.
To determine the average inhibitory concentration (IC 50 ) for each molecule, the ligand efficiency equation vide infra was used: where N is the number of non-hydrogen atoms, and pIC 50 is the negative logarithm of IC 50 . Knowing the theoretical LE value through previous calculations and the molecular docking ∆G values, the formula was adjusted to obtain the pIC 50 value and later the IC 50 value by the following equation: pIC 50 = LE × N/1.4 (4)

Conclusions
With this study, the probable region of the active theoretical site of the crystal structure of Spike-ACE2 was determined. Using hydroxychloroquine as a standard, we obtained the pharmacophoric model from antiparasitics. Tanimoto analysis and hierarchical analysis studies raised the hypothesis that structures with the same chemical characteristics would be great candidates for SARS-CoV-2 inhibitors.
The pharmacophores generated a list of chemical structures obtained by virtual screening that were evaluated. The molecules approved in the ADME/Tox tests showed excellent synthetic viability values, concluding that they are easy to obtain through synthesis. The pharmacokinetic and pharmacological properties showed better parameters than the control drug, and none showed carcinogenic potential, performing even better than hydroxychloroquine.
All studied molecules showed better results than hydroxychloroquine in molecular docking tests, with ligand_003 having a strong binding affinity of −8.645 kcal·mol −1 . The lipophilicity and solubility data showed favorable values for oral administration. The theoretical IC 50 values showed that these molecules are promising for in vitro studies and the chemical descriptors demonstrate great stability for the molecules.