Viroinformatics-Based Analysis of SARS-CoV-2 Core Proteins for Potential Therapeutic Targets

SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) is a novel coronavirus for which no known effective antiviral drugs are available. In the present study, to accelerate the discovery of potential drug candidates, bioinformatics-based in silico drug discovery approaches are utilized. We performed multiple sequence alignments of the Spike (S) protein with 75 sequences of different viruses from the Orthocoronavirinae subfamily. This provided us with insights into the evolutionarily conserved domains that can be targeted using drugs or specific antibodies. Further, we analyzed the mechanism of SARS-CoV-2 core proteins, i.e., S and RdRp (RNA-dependent RNA polymerase), to elucidate how the virus infection can utilize hemoglobin to decrease the blood oxygen level. Moreover, after a comprehensive literature survey, more than 60 antiviral drugs were chosen. The candidate drugs were then ranked based on their potential to interact with the Spike and RdRp proteins of SARS-CoV-2. The present multidimensional study further advances our understanding of the novel viral molecular targets and potential of computational approaches for therapeutic assessments. The present study can be a steppingstone in the selection of potential drug candidates to be used either as a treatment or as a reference point when designing a new drug/antibody/inhibitory peptide/vaccine against SARS-CoV-2.


Introduction
In the past two decades, the world has faced several infectious disease outbreaks, such as influenza A (H1N1), SARS, MERS, Ebola, and Zika virus. These have had an and similar pathogens. The present report provides insights into the structural motifs of the SARS-CoV-2 core proteins and assessed therapeutic potency of existing drug candidates against SARS-CoV-2.

Investigation of Conservative Motifs of SARS-CoV-2 Spike Protein and Phylogenetic Analysis
To investigate the SARS-CoV-2 Spike protein's conservative motifs and submit it to phylogenetic analysis, we used the Blast Molecular Evolutionary Genetics Analysis (MEGA 10.0.5) tool. In search of unique conservative regions in the Spike protein throughout its evolution, we performed alignment of the amino acid sequence of the Spike protein with the BLASTP database. We then selected 75 subject sequences of different viruses from the Orthocoronavirinae subfamily, which had the highest BLASTP score. The maximum number of possible alignments in one run determined the number of sequences. Later, we ran the alignment of these 75 sequences from the database using the MEGA-X software's MUSCLE program [19]. Further, a phylogenetic tree was generated by using the maximum likelihood method and the JTT matrix-based model [20]. We then selected the tree with the highest log-likelihood score for the top 17 virus strains that are phylogenetically closely related to the novel coronavirus. The percentage of trees in which the associated taxa clustered together is shown next to the branches that were obtained. This analysis involved 17 amino acid sequences. There were a total of 1728 positions in the final dataset. The phylogeny was tested using the resampling bootstrap methodology with a bootstrap value of 100 [19].
The five selected species that underwent detailed alignment evaluations are as follows:

In Silico Protein-Protein Interaction
We performed the protein-protein docking using ClusPro (v2.0) server [21]. We submitted macromolecules in .pdb format to the ClusPro (v2.0) server and downloaded the most favorable interaction model for further analysis in Discovery Studio Visualizer (DSV) (v2019 client; Dassault Systèmes BIOVIA). We docked the SARS-CoV-2 Spike protein with both the tetramer and monomer molecules of hemoglobin (1A3N.pdb) and the heme group (extracted from 1A3N.pdb), in ClusPro server and AutoDockTools (v1.5.6) software, respectively. LigPlot + v.2.2 software was used to generate a 2D plot of interrace amino acids at the Hemoglobin-Spike interaction pocket [22].

Finding the Most Interacting Motifs of Spike Protein with Antiviral Drugs and Heat-Map Representation of Motifs
We manually investigated the binding pocket of Spike protein and RdRp with the docked drug ligands. We selected antiviral FDA approved drugs, flavonoids, and antibiotics, which are readily available on the market. In addition to that, these drugs are under investigation for their potency in the treatment of COVID-19. We selected the drugs from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), which is a national library of medicines managed by National Institute of Health, USA. We divided the whole protein sequence of Spike and RdRp proteins in, respectively, 12 and 10 equally spanned motifs. Based on the number of interactions with the docked drug ligand database, we ranked each motif and further normalized the number of interactions in each motif to compare their potential as a target for the therapeutic drugs and their significance for future drug development against SARS-CoV-2. In continuation, based on the normalization of the total number of interactions for each motif, heat-maps were produced using DSV. Similarly, we investigated the frequency and location of specific amino acids of Spike and RdRp in each motif, which are important for the interaction with various drug ligands.

Docking of Antivirals, Antibiotics, Antiparasitics, Flavonoids and Vitamins with Spike and RdRp Protein of SARS-CoV-2
Three-dimensional models of SARS-CoV-2 Spike (QHD43416.pdb) and RdRp (QHD43 415_11.pdb) proteins were downloaded from Zhang Lab I-TESSAR online library in .pdb format. We downloaded 3D structures of drugs and small molecules for docking from PubChem in .pdb format [23]. Literature suggested, there is no significant change in the structure of these drug molecules after getting activated inside the liver, which makes them a suitable candidate for the docking study. Docking was performed with the AutoDockTool (v1.5.6) software [24]. There was no change observed in these drug molecules' structure after being activated inside the liver, which makes them a suitable candidate for the docking study. Docking was performed with the AutoDockTool (v1.5.6) software [24]. All models were converted to "pdbqt" format after deleting water molecules and computing the Gasteiger algorithms for the macromolecules, whereas, in ligands, we choose torsion angles ≤ 9. For developing the grid box for SARS-CoV-2 Spike protein, we set the following metrics: X-dimension: 104, Y-dimension: 116, Z-dimension: 126, Spacing: 1.000, Offset X: 2.917, Offset Y: 1.450 and Offset Z: 0.194 and grid was set using Autogrid4. Docking was performed with Autodock4; for each docking, we set our macromolecules' Rigid filename and docked with Genetic algorithm and output set to Lamarckian GA (4.2). After the docking simulation, 10 docking models were analyzed and ranked by their energy binding ability (∆G). H-bonds were also built in the docking models, and the best model of each docking was saved as a ".pdbqt" format file. Docking models were further analyzed for the discovery of a binding pocket in Discovery Studio Visualizer (DSV) (v2019 client; Dassault Systèmes BIOVIA) software. Interactions were built by selecting all favorable interaction types intermolecularly. Furthermore, 2D diagrams were visualized through DSV.

Phylogenetic Analysis of SARS-CoV-2 Spike Protein and Retrieval of Therapeutically Important Conservative Motifs of Evolutionary Significance
We ran the sequence alignment of 75 subject sequences from the database using MEGA-X software and developed a phylogenetic tree through BLASTP. The subject sequences suggested by BLASTP are related to the Spike (S) proteins from MERS, murine hepatitis virus (MCoV), human coronavirus (HCoV), feline infectious peritonitis virus (FIP-CoV) and SARS-CoV-1 ( Figure 1). The initial results suggest that the SARS-CoV-2 Spike protein shares similarities with SARS-CoV-1, whereas a common ancestor connects MERS and SARS-CoV-2. The majority of similar sequences aligned with SARS-CoV-2 Spike protein are from the Betacoronavirus of the subspecies of Sarbecoviruses, Merbecoviruses, and Embecoviruses; conversely, the minority of them come from Alphacoronaviruses.
Similarly, we identified conserved regions of Spike protein through the sequence alignment among the subjected 75 sequences. Starting from S1 N-terminus Domain (NTD), 10 out of 75 subjected sequences match with its first portion at the 1-200 amino acid (aa) range, whereas 24 out of 75 sequences match with the S1-NTD's last portion at 201-300 aa. A moderate amount of similarity was observed at the S1 C-terminus Domain (CTD) side with 33 similar sequences at around 300-510 ( Figure 2). Moreover, the highest conservation among subjected sequences was observed in the center of the S2 portion of the protein at around 900-1000 aa in 42 out of 75. To strengthen these results, we aligned five Spike proteins among several species related to coronavirus and investigate their conservation in five core domains (S1-NTD, S1-CTD, PCS, S2-HR1, and S2-HR2). The five additional coronavirus species sequences that were subjected to detailed alignment were chosen by their displayed values in the phylogenic tree. As a reference value, we chose the spike protein of novel coronavirus (PDB ID: 6VXX), then we chose three species possessing the highest values and two species possessing the lowest values. This will help us to include a spectrum of proteins that relates to an ancestor protein. The results suggest that the domains of S1-NTD, HR1, and HR2 are highly conserved ( Figure 2). Similarly, we identified conserved regions of Spike protein through the sequence alignment among the subjected 75 sequences. Starting from S1 N-terminus Domain (NTD), 10 out of 75 subjected sequences match with its first portion at the 1-200 amino acid (aa) range, whereas 24 out of 75 sequences match with the S1-NTD's last portion at 201-300 aa. A moderate amount of similarity was observed at the S1 C-terminus Domain (CTD) side with 33 similar sequences at around 300-510 ( Figure 2). Moreover, the highest conservation among subjected sequences was observed in the center of the S2 portion of the protein at around 900-1000 aa in 42 out of 75. To strengthen these results, we aligned five Spike proteins among several species related to coronavirus and investigate their conservation in five core domains (S1-NTD, S1-CTD, PCS, S2-HR1, and S2-HR2). The five additional coronavirus species sequences that were subjected to detailed alignment were chosen by their displayed values in the phylogenic tree. As a reference value, we chose Furthermore, we could identify several conserved regions spanning throughout the Spike protein among the selected viral strains. There are several "clusters" of conserved regions spanning from the beginning S1-CTD until its end at the 330 to 529 aa range. This means that S1-CTD can be divided into seven clusters of conserved domains, composed of around 25 amino acids each, according to the MEGA-X software alignment ( Figure S1). This expands our understanding of nature of the Spike protein Binding Domain to ACE2 and other receptors. Moreover, two very large evolutionarily conserved domains lie at the center of the S2 region within the amino acid ranges of 908 to 1003 and 1163 to 1211, indicating that the S2 part of Spike protein seems to be heavily conserved in comparison to its S1 counterpart. Taken together, the S1-CTD and the final portion of the S2 protein are relatively conserved throughout the species and both serve as effective references for the universal drug discovery against most of the strains of SARS-CoV-2 ( Figure 2). Each lane consists of 6 sequences subjected to the alignment (n). Conservation exists in higher levels at S2-HR1 and HR2, moderate conservation occurs at S1-NTD and PCS, whereas at S1-CTD there are many hyper-variable regions.

Figure 3. (A)
Three-dimensional representation of the interaction between the Spike protein (in yellow color) and hemoglobin monomer (in magenta color) and detailed representation of the binding area. (B) Two-dimensional interaction plot was formed via LigPlot. Interaction may occur at S1-NTD with the referred amino acids from each side. In the green dashed line, hydrogen bonds (B) Two-dimensional interaction plot was formed via LigPlot. Interaction may occur at S1-NTD with the referred amino acids from each side. In the green dashed line, hydrogen bonds are represented and radian arches represent weaker electrostatic bonds that participating in the interaction.
Interestingly, the ARG-237 hydrogen bond persists in both simulations (monomer and tetramer), thus indicating a possible binding area that is more favorable for the interaction in terms of binding efficiency and spatial pairing. Moreover, the region surrounding the above-mentioned amino acid might facilitate the binding of the heme group to the S1-NTD region of the Spike protein. Simulation of this interaction yielded a binding energy score of −5.4 kcal/Mol. The results of the interaction show several Van der Waals bonds surrounding the heme group, whereas there are six Pi-(Anion, Sigma, -Pi Stacked and Alkyl) bonds and one hydrogen bond bridging with the Spike protein ( Figure 4A). Amino acid residues TYR-38, LYS-41, PHE-43, LYS-206, GLU-224, PRO-225, LEU-226, VAL-227, ASP-228, and THR-284 of Spike protein interact with the heme group directly. The core amino acid of this interaction is GLU-224, which might possess the ability to directly bind the heme iron with a Pi-Anion bond, thus preventing free oxygen-binding in this area ( Figure 4B). bonds surrounding the heme group, whereas there are six Pi-(Anion, Sigma, -Pi Stacked and Alkyl) bonds and one hydrogen bond bridging with the Spike protein ( Figure 4A). Amino acid residues TYR-38, LYS-41, PHE-43, LYS-206, GLU-224, PRO-225, LEU-226, VAL-227, ASP-228, and THR-284 of Spike protein interact with the heme group directly. The core amino acid of this interaction is GLU-224, which might possess the ability to directly bind the heme iron with a Pi-Anion bond, thus preventing free oxygen-binding in this area ( Figure 4B).

Most Important Conservative Motifs of SARS-CoV-2 Spike Protein as a Target for the Development of Therapeutic Drugs
Most of the antiviral drugs interact with Spike S1-NTD, especially between the residues 1-100 and 201-300 aa ( Figure 5A) and with most common interacting residue locations at VAL-3, PHE-4, VAL-6, PHE-58, and PRO-82. On the other hand, antiviral interaction occurs with Spike S2, lying between the residues 801-1000 and 1201-1300 aa. In particular, the amino acids ASN-824, VAL-826, THR-827, LEU-945, LYS-1205, TYR-1209, PRO-1213, and TRP-1217 are very important due to their participance in most of the interactions (Table S1). Simulations suggest a preference of antiviral drugs to bind around the S1-NTD region and S2 (A); similar results were observed in flavonoid interactions (D). Docking of antibiotic drugs exhibits a preference in the S1-NTD region only (B). Antiparasitic drugs interestingly showed a higher binding frequency in S2 of the protein (C). From the perspective of the vitamins, they were shown to bind in both the S1-NTD/-CTD region and S2 (E). Simulations suggest a preference of antiviral drugs to bind around the S1-NTD region and S2 (A); similar results were observed in flavonoid interactions (D). Docking of antibiotic drugs exhibits a preference in the S1-NTD region only (B). Antiparasitic drugs interestingly showed a higher binding frequency in S2 of the protein (C). From the perspective of the vitamins, they were shown to bind in both the S1-NTD/-CTD region and S2 (E).
Furthermore, in the case of antiparasitic drugs, we found that the residues of S2 Spike between 901-1000 and 1101-1200 aa are very important, especially the locations LYS-1205 and TYR-1209, which are similar to antivirals and antibiotics ( Figure 5C). However, we also found interaction at S1-NTD between the residues 1-101 and 200 aa, which suggests the therapeutic importance related to S1 Spike protein. Additionally, a study with flavonoids revealed that interactions mainly happen at S1-NTD between the residues 1-100, especially at the locations VAL-6 and ARG-21, and at S2 between the residues 701-900, such as LYS-790, LYS-811, ALA-893 and LEU-1224 ( Figure 5D). Similarly, the interaction with vitamins C and D lies between the residues 1-100, 301-400, and 501-600 amino acids of Spike S1, most importantly VAL-47, PHE-318, and THR-630, and between the 701-1000 residues of Spike S2, especially the residues TYR-741, ILE-742, CYS-743, and GLY-744 ( Figure 5E; Table S3).

Most Important Motifs of RdRp Protein as a Target for Therapeutic Drugs
Most of the drugs targeting RdRp mainly interact between the amino acid residues 101-200 and 701-800, especially THR-51, ASN-122, PHE-133, LYS-182, CYS-760, THR-761, ALA-771, GLU-780 ( Figure 6A,B). Additionally, we also found a subtle interaction between 1-100 and 801-900 amino acid residues ( Figure 6B; Table S4).   Table S4. Binding energies are organized here (A) and are arranged from weakest to strongest based on the molecule they are interacting with. These lists of drugs and their binding energies are separated along the central axis based on the viral protein they are interacting with, namely RdRp on the left and the Spike protein on the right-hand side. Additionally, quantitative analysis and graphical representation of the locations at which the molecules are interacting with the SARS-CoV-2 RdRp (B) can be observed.

Discussion
In the search for unique conservative regions in the SARS-CoV-2 Spike protein, we performed an alignment of the amino acid sequence of the Spike protein with the BLASTP database ( Figure 1; Figure S1). Spike S1 contains a receptor binding domain (RBD) (AA: 449-510) actively binding with ACE2 receptors (Figure S2), which contains several residues evolutionary conserved and dedicated for the interaction between ACE2 and SARS-CoV strains including SARS-Cov-2 [25][26][27]. In the current study, alignment showed that conservation exists in higher levels at S2-HR1 and HR2, moderate conservation at S1-NTD, and PCS, whereas at S1-CTD there are many hypervariable regions (Figure 2). This might suggest that drugs which tend to bind to those regions are more effective in addressing a bigger spectrum of coronaviruses. Interestingly, most of the antiviral top hit drugs, such as Indinavir (∆G = −9.8 kcal/Mol) and Nelfinavir (∆G = −9 kcal/Mol), bind near the center of S2 (Figure 2; Table S1). A similar pattern was observed with antiparasitic drugs Ivermectin B1a/b (∆G = −9.16/−8.86 kcal/Mol) and vitamin D (∆G = −5.52 kcal/Mol), whereas, on the other hand, top hit antibiotic drugs did not show preferences at the S2 binding area of the Spike protein (Tables S2 and S3). From this perspective, the above-mentioned drugs could possess the ability to physically inhibit the assembly of the Spike protein as well as the viral entry in a pan-CoV targeted manner.
Shuai and colleagues showed that targeting the HR1 domain of Human-CoV Spike protein, which is in the area of S2, showed a significantly reduced viral entry in mouse models [28]. Our docked drugs are also bind in the same region of Spike, supporting the results of earlier reports [28]. Targeting this specific domain with a small molecule could be a pan-CoV target, as its conservation is high among coronavirus species. Turning the attention to conserved areas of Spike protein and not to the hypervariable S1-CTD will lead to a drug that not only has great inhibitory potential but also will possess a broad target activity against multiple coronavirus infection.
Furthermore, taking into consideration another spotlight region of the Spike protein, the polybasic furin cleavage site (PCS), which is located in the boundaries of S1 and S2, was under observation [1]. Throughout our docked library, only one molecule could bind to that area. Vancomycin, which exhibited the lowest binding energy (∆G = −10.2 kcal/Mol) from all our docked drug categories, is the only drug that binds within the region of PCS in an approximate position of amino acid range between 654 and 692. The functional activity is still unknown for SARS-CoV-2 PCS; pathogenicity and transmissibility need to be elucidated as well. Emerging research referred to PCS as an enhancer of cell-cell fusion, whereas cleavage of it results in intra-species infection, with a great example being MERS-CoV [29]. Thus, it is very important to evaluate drugs that can bind to that area of the Spike protein since it might result in desirable outcomes supporting the fight against SARS-CoV-2.
Recently, many anecdotal clinical trial reports have suggested that patients affected by COVID-19 not only have impaired respiratory function but also low oxygen count, causing respiratory distress [7,15,16]. This phenomenon led us to test for a possible binding event between the SARS-CoV-2 Spike protein and the heme group of the hemoglobin molecule. That binding might greatly reduce the oxygen affinity to hemoglobin and increase its degradation rate, as many COVID-19 survivors showed a much lower hemoglobin count compared with those not infected [17,18]. It was also revealed that alternatives currently being considered for COVID-19 treatment by increasing hemoglobin production and increasing hemoglobin availability for oxygen binding and by causing hyperventilation with associated increasing levels of oxygen and decreasing levels of carbon dioxide in the blood may significantly ameliorate COVID-19 respiratory symptoms [30,31].
In line with our findings, Spike protein has the potential to interact with the tetramer and monomer of hemoglobin (Figure 3; Figure S3). Here, amino acid GLU224 might be the key target as it plays a crucial role in binding directly with the iron particle of heme group ( Figure S3). This is the first study which shows that the Spike protein of SARS-CoV-2 might interact with not only the hemoglobin molecule but with the heme group as well. These results will shed light on the, yet unknown, intrinsic mechanism of SARS-CoV-2 Spike protein and the low oxygen count of patients and survivors, while, at the same time, while serving as will serve as a basis for the development of new drugs.
Based on the simulated interaction of SARS-CoV-2 Spike protein with crucial antiviral, antibacterial, antiparasitic, flavonoids and vitamins, we found a few striking readouts. Initially, we focused on 47 different antivirals already approved and broadly available in the clinical practice suggesting these have the potential of a viral antagonistic effect. Given the fact that the Spike protein is evolutionarily conserved, it is a prime location for interaction and inhibition of the virus in a manner that is ubiquitous across all potential strains. Our results show a wide range of effectiveness of antivirals drugs, from Saquinavir at the lowest (∆G = −0.23 kcal/Mol) to Indinavir at the highest (∆G = −9.8 kcal/Mol) ( Figure S4). Drugs with interaction energy greater than or equal to the Japan-made Favipiravir (∆G = −5.2 kcal/Mol) might be prime candidates to be potential treatments, as it has already been observed to have mild inhibiting effects on COVID-19 patients [32]. Further examination of the results revealed that the drug molecules interacted with the Spike protein monomer at very specific locations ( Figure 5A). Most notable are the β-sheets within the 1-100 residues range having over 150 different molecular interactions. This is followed by the α-helix within the 801-900 aa residues range yielding over 100 unique interactions. Noting the location of the 1-100 aa range, drugs that interact with that binding pocket will most likely not have a large-scale effect on the protein's functionality, as it only serves as the N-terminal domain. Conversely, the 801-900 aa region is located at the point on the Spike monomer that is key to the formation of the Spike trimer; binding to this region would impair SARS-CoV-2's ability to form the trimeric structure and overall inhibit the production of viable virions [11,12].
Researchers have shown that the combined use of azithromycin and hydroxychloroquine suppressed the growth of SARS-Cov-2 virus in vitro [33,34]. Therefore, similar to the examination of antiviral drugs, we investigated the binding properties of 21 prominent antibiotics to the SARS-CoV-2 Spike protein. Unfortunately, the majority of the molecular interactions observed occur in the N-terminal domain of the Spike protein ( Figure 5B), which indicates that they might yield little to no result in disrupting the virus. In conclusion, antibiotic treatment may not play an important role in combinational therapies.
The prospective anti-parasitic, Ivermectin and hydroxychloroquine, the two FDAapproved drugs that were previously shown to have broad-spectrum antiviral activity [33,35]. Leon Caley et al. showed that a single addition of Ivermectin to Vero-hSLAM cells 2 h post-infection with SARS-CoV-2 significantly reduced viral RNA [36]. Ivermectin, therefore, warrants further investigation for possible benefits in humans. Similarly, in our study, we found that Ivermectin B1a and B1b have low energy interactions (∆G = −9.16 kcal/Mol and ∆G = −8.86 kcal/Mol respectively). Additionally, these interactions were observed to occur primarily on the α-helices existing in the 901-1200 aa region of the Spike monomer ( Figure 5C), given the fact that this region is involved in the formation of the Spike trimer, specifically the base of Spike where protein is integrated into the virion, which is conserved among coronavirus species [2,11,12]. Similarly, the 4-aminoquinoline antimalarials chloroquine and hydroxychloroquine have been investigated (Gautret, Lagier et al., 2020), which further support the interaction we found between hydroxychloroquine and the SARS-CoV-2 Spike protein (∆G = −2.85 kcal/Mol) ( Figure S4).
In this line, flavonoids were shown to play a crucial role against different emerging viruses. Tetramethoxyflavone from elderberry extract inhibited human influenza A (H1N1) infection in vitro according to a study of Roschek and Fink [37]. Moreover, flavonoids were utilized against the infectious bronchitis virus (IBV), a pathogenic chicken coronavirus [38]. Chen and colleagues tested non-cytotoxic, crude ethanol extracts of Sambucus nigra fruit for anti-IBV activity and showed inhibition of the virus by compromising their envelopes. They also suggest that future studies using S. nigra extract to treat or prevent IBV or other coronaviruses are warranted [39]. It has been suggested that flavonoids such as herbacetin, rhoifolin, and pectolinarin were found to effectively block the enzymatic activity of SARS [40,41]. In our study, we focused on five different flavonoids ( Figure S4) that were proven to have some effect with other coronaviruses. Although flavonoids did not exhibit high ∆G values, they showed unique binding areas. From the binding areas, we could roughly identify the preferred binding regions on the Spike protein ( Figure 5). Total interactions per segment results suggest a tendency of flavonoids to bind in both S1-NTD and S2. Even though most of the screened drugs and small molecules bind with the S1-NTD region, flavonoids showed an increased number of interactions in the conserved S2 region too. This suggests that they might possess the ability not only for structurally inhibiting the assembly of the Spike trimer but also for binding into the lumen of the trimeric protein and thus inhibiting the viral fusion with the host cell in a variety of coronaviruses.
Earlier, researchers claimed that vitamin C and vitamin D have a beneficial effect on lower respiratory tract diseases [42,43]. In the light of the evidence from various groups regarding the beneficiary effect of these vitamins, we investigated their interaction with the Spike protein. There are two forms of vitamin D-cholecalciferol (vitamin D3) and ergocalciferol (vitamin D2) [44]. Here, we found that vitamin C (∆G = −2.95 kcal/Mol) and vitamin D3 (∆G = −5.52 kcal/Mol) can directly interact with the Spike protein of SARS-CoV-2 (Table S3), which further reaffirmed the possibility of their use as a drug for the treatment of COVID-19. Recently, researchers reported an association between reduced levels of vitamin D and mortality cases caused by COVID-19 in various countries. These results suggest that a higher mean level of vitamin D could results in lower mortality rates [45].
Inhibition of the RNA-dependent RNA polymerase (RdRp) employed by RNA viruses during their replication process [46,47] has been an effective treatment strategy. Japanese anti-influenza drug Favipiravir (sold as Avigan) showed promise in the treatment of COVID-19 patients. However, as the pandemic evolved, it became clear that this drug was only effective if administered at the early stages of coronavirus infection [46]. This suggests that the inhibition of RdRp works primarily as a preventative measure rather than a primary form of treatment against RNA-based viruses. Further, our findings suggest the efficacy of Beclabuvir, in particular, to be used against SARS-CoV-2. Not only does it interact significantly with the viral RdRp, but its binding energy to the Spike protein exceeds that of a significant number of the drugs tested in this study. This suggests that designing a drug which targets similar residues would have a stronger inhibitory effect on the virus than that of Favipiravir and provide additional protection in the form of disrupting the Spike protein's function. Using Favipiravir as a standard at ∆G = −3.09 kcal/Mol for the ability to exhibit a mild clinical success, three RdRp-targeting drugs are immediately identifiable as potentially more successful: Ribavirin (∆G = −4.20 kcal/Mol), Galidesivir (∆G = −4.38 kcal/Mol), and Beclabuvir (∆G = −5.63 kcal/Mol). Additional examination of where these molecular interactions take place revealed that they are primarily concentrated in the 101-200 and 701-800 aa regions ( Figure 6B). Similar to the Spike protein, the 701-800 aa regions in particular consist of a series of α-helices that are essential to the RdRp's functionality; blocking them would, therefore, inhibit the use of RdRp and essentially shut down the virus's ability to replicate its genome.
Taken together, analysis of the most interacting motifs of Spike and RdRp, along with evolutionary conserved motifs of SARS-CoV-2 Spike, provided advanced understanding of the computer-aided drug and antibody/vaccine designing for the treatment of COVID-19. Our analysis suggests there is a high possibility of mutation in the S1 region of the Spike protein; in contrast, S2 remains more conservative during the evolution. Therefore, the use of a combination of conserved peptide sequences of Spike protein could be a strategy to develop an effective vaccine. In addition, designing the new drugs selectively targeting these conservative motifs might help to find out a potential cure for the COVID-19. However, prior in vitro/in vivo investigation is warranted to ensure the clinical application of the results of the present study. Conclusively, designing inhibitory peptides blocking the receptor-binding domain of Spike protein interacting with hemoglobin/heme that we reported in the current study could help to control the hypoxia condition during the infection of SARS-CoV-2, which might significantly improve the survival of patients suffering from COVID-19.

Supplementary Materials:
The following are available online at https://www.mdpi.com/2073-4 468/10/1/3/s1. Table S1: List of the ten top hits of antiviral drugs against SARS-CoV-2 Spike protein by order of lowest to highest binding energy value. Indinavir seems to exhibit the lowest binding energy (−9.8 kcal/Mol) when docked with Spike protein, whereas Cobicistat exhibits the highest binding energy (−6.37 kcal/Mol). Table S2: List of the ten top hits of antibiotics drugs against SARS-CoV-2 Spike protein by order of lowest to highest binding energy value. Vancomycin seems to exhibit the lowest binding energy (−10.2 kcal/Mol) when docked with Spike protein, whereas Levofloxacin exhibits the highest binding energy (−5.11 kcal/Mol). Table S3: List of the top hits antiparasitic drugs, flavonoids, and vitamins against SARS-CoV-2 Spike protein by order of lowest to highest binding energy value. Ivermectin B1a seems to exhibit the lowest binding energy (−9.16 kcal/Mol) when docked with Spike protein, whereas hydroxychloroquine exhibits the highest binding energy (−2.85 kcal/Mol). In terms of flavonoids, tetramethoxyflavone exhibits the lowest binding energy (−4.91 kcal/Mol) and Gallocatechin exhibits the highest binding energy (−2.88). Moreover, vitamin D showed the lowest binding energy (−5.52 kcal/Mol) and vitamin C the highest binding energy (−2.95 kcal/Mol). Table S4: A comparative analysis between SARS-CoV-2 RdRP binding affinities and spike protein binding affinities in antivirals which primarily target viral RdRP. These comparisons include binding energy, the specific amino acids that are interacting with the drug, and the types of bonds formed during this interaction. Figure S1: The tree was constructed based on neighbor-joining and BIONJ algorithms using MEGA 6.0. The topological stability of the NJ tree was achieved by running 100 bootstrapping replications. Bootstrap values (%) are indicated by numbers at the nodes and the length of the branch (represents the evolutionary time between two nodes). This analysis involved 75 amino acid sequences. Figure S2: Interaction of the trimeric SARS-CoV-19 Spike protein with the dimeric ACE2 receptor using the ClasPro server. Interaction is happening at the S1-CTD of the Spike protein (A). Interacting amino acids of each protein responsible for the interaction visualized by Ligplot. Green dashed lines between the opposing amino acids represent hydrogen bonds, whereas arches represent other forces, i.e., Van der Waals forces (B). Figure S3: (A) Three-dimensional representation of the interaction between the Spike protein (in yellow color) and hemoglobin tetramer (in magenta color). (B) Two-dimensional interaction plot was formed via LigPlot. Interaction may occur at S1-NTD with the referred amino acids from each side. The green dashed lines represent hydrogen bonds between the interactive amino acids. Figure S4: Circular tree represents the order of antiviral drugs (A) and antibiotics, antiparasitic, flavonoids, and extra ligands (B) based on the binding energy (∆G in kcal/Mol).
Author Contributions: L.A. and T.P. developed the idea, designed the study, performed the simulations, data analysis, and contributed equally to the manuscript. S.E., C.M., R.K.B., M.-H.C. and R.S.K. contributed to the development of the idea, and contributed to writing and editing the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding:
The study was financially supported by Ministry of Science and Technology (Grant No. MOST 108-2218-E-033-005-MY2). Authors utilized no grant support for the current research. Only the publication fee will be supported by the grant.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article [and/or] its supplementary materials. Generated raw data supporting the findings of this study are available from the corresponding author [LA] on request.