Antiviral Efficacy of Selected Natural Phytochemicals against SARS-CoV-2 Spike Glycoprotein Using Structure-Based Drug Designing

SARS-CoV-2 is a highly virulent coronavirus that first surfaced in late 2019 and has since created a pandemic of the acute respiratory sickness known as “coronavirus disease 2019” (COVID-19), posing a threat to human health and public safety. S-RBD is a coronaviral protein that is essential for a coronavirus (CoV) to bind and penetrate into host cells. As a result, it has become a popular pharmacological target. The goal of this study was to find potential candidates for anti-coronavirus disease 2019 (COVID-19) drugs by targeting severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) S-RBD with novel bioactive compounds and molecular interaction studies of 15,000 phytochemicals belonging to different flavonoid subgroups. A spike protein crystal structure attached to the ACE2 structure was obtained from the PDB database. A library of 15,000 phytochemicals was made by collecting compounds from different databases, such as the Zinc-database, PubChem-database, and MPD3-database. This library was docked against a receptor binding domain of a spike glycoprotein through the Molecular Operating Environment (MOE). The top drug candidates Phylloflavan, Milk thistle, Ilexin B and Isosilybin B, after virtual screening, were selected on the basis of the least binding score. Phylloflavan ranked as the top compound because of its least binding affinity score of −14.09 kcal/mol. In silico studies showed that all those compounds showed good activity and could be used as an immunological response with no bioavailability issues. Absorption, distribution, metabolism, excretion and a toxicological analysis were conducted through SwissADME. Stability and effectiveness of the docked complexes were elucidated by performing the 100 ns molecular dynamic simulation through the Desmond package.


Introduction
Coronavirus (CoV) is a virus encircled by a single-stranded, positive-sense RNA genome that is believed to produce human respiratory infections. The world has been grasped in a pandemic by the attack of the acute respiratory disease coronavirus (SARS-CoV-2) that was established in 2019. The world's culture, economy, healthcare and infrastructure of its population has been challenged by this pandemic [1,2]. In December 2019, an unusual novel and viral pneumonia originated in Wuhan, China. It has been recognized as a zoonotic disease similar to SARS-CoV and MERS-CoV and, hence, given the name nCoV-2019. In January 2020, a patient affected by the coronavirus was observed using the throat swab test by the Chinese center for disease control and prevention and, consequently, it was identified as nCoV-2019 by the World Health Organization [3,4]. This epidemic illness caused by SARS-CoV-2 was announced as coronavirus disease   [4]. Once a person has become infected with corona virus, it takes 2-14 days for the appearance of symptoms [5]. A period of 5-6 days is roughly considered as the average incubation receptor-binding domain, we first studied the spike receptor-binding domain and the H-ACE2 complex interfaces. Some potent spike receptor-binding domain residues significantly bonded to the ACE2, leading to the stability of complexity. Figure 1. 3D-Visulization of the spike glycoprotein receptor with human ACE2 enzyme.

Database Screening and Docking Study
The spike protein was docked against a ligand molecular library of 15,000 phytochemicals. Docking was a reliable procedure which showed the relationship between phytochemicals and the SARS-CoV-2 spike protein receptor-binding domain. MOE software was used to perform molecular docking in order to find compounds with the best binding residue-to-receptor interaction. The RMSD value, the residue binding with ligands and the docking score were used to evaluate compounds. The top six ranked docking poses were picked from 15,000 docked molecules. Compounds with a low score, such as rmsd3, and residues with the most interactions were chosen. The lowest binding energy of these selected phytochemicals ranged from −14.09 to −12.19 (kcal /mol). In order to obtain a precise sense of the receptor-ligand interaction with the highest docked complexes, the minimum binding energy and scoring function of each docked complex were shown in Table  1.  Figure  3.

Database Screening and Docking Study
The spike protein was docked against a ligand molecular library of 15,000 phytochemicals. Docking was a reliable procedure which showed the relationship between phytochemicals and the SARS-CoV-2 spike protein receptor-binding domain. MOE software was used to perform molecular docking in order to find compounds with the best binding residue-toreceptor interaction. The RMSD value, the residue binding with ligands and the docking score were used to evaluate compounds. The top six ranked docking poses were picked from 15,000 docked molecules. Compounds with a low score, such as rmsd3, and residues with the most interactions were chosen. The lowest binding energy of these selected phytochemicals ranged from −14.09 to −12.19 (kcal /mol). In order to obtain a precise sense of the receptorligand interaction with the highest docked complexes, the minimum binding energy and scoring function of each docked complex were shown in Table 1. Two-dimensional maps of these interactions were analyzed using the MOE LigX tool. ChimeraX was used to display docked complexes. The LigX interaction diagram showed that Phylloflavan and Milk thistle were found to bind with the spike receptor binding domain protein with a score of −14.09 kcal/mol and −13.

Drug-Likeness/ADMET Profiling
To analyze the drug-likeness qualities of the top compounds, a drug scanning was performed using the Molinspiration service. Lipinski's Rule of Five has become a norm; this rule depicts important drug features such as pharmacokinetics, interactions and metabolism in the human body, as well as their excretion. The selected compounds had no violations of Lipinski's five rules and had significant drug-like features, such as molecular weight (Table 2).

Drug-Likeness/ADMET Profiling
To analyze the drug-likeness qualities of the top compounds, a drug scanning was performed using the Molinspiration service. Lipinski's Rule of Five has become a norm; this rule depicts important drug features such as pharmacokinetics, interactions and metabolism in the human body, as well as their excretion. The selected compounds had no violations of Lipinski's five rules and had significant drug-like features, such as molecular weight ( Table 2).

Drug-Likeness/ADMET Profiling
To analyze the drug-likeness qualities of the top compounds, a drug scanning was performed using the Molinspiration service. Lipinski's Rule of Five has become a norm; this rule depicts important drug features such as pharmacokinetics, interactions and metabolism in the human body, as well as their excretion. The selected compounds had no violations of Lipinski's five rules and had significant drug-like features, such as molecular weight (Table 2).   ADME and AdmetSAR were used to assess a number of pharmacokinetic variables. The ADME and toxicity of the top therapeutic candidate medications can be estimated using pharmacokinetic characteristics. Table 3 shows the ADMET properties of derived phytochemicals for both targets. Due to poor pharmacokinetic qualities and toxicity, many medicines do not exploit this mechanism in their development. Early drug discovery relies on high-performance and quick ADMET profiling investigations to identify active lead compounds. SwissADME was also utilized to test the ADMET profiles of prospective compounds to validate the drug likeliness (Table 3).

Biochemical Classification of Idetified Compounds
All the compounds' classifications and their use against diseases were previously explained in Table 4. All structures and chemical properties, including taxonomy and its biological classification, were explained for a better understanding of the efficacy of our predicted drug candidates.

Energy Calculations
The binding free energy (deltaG bind) of inhibitory drugs with a strong potential to inhibit the activity of the SARS-CoV-2 spike protein was calculated using MMGBSA. To perform the MMGBSA energy estimates, docking complexes with high energy function scores were obtained. By minimizing protein-ligand complexes, salvation energy and surface area energy, total free binding energy was determined. The results revealed that the compound name with ID had a high level of stability, with a total binding free energy of −30.35 kcal/mol for the Phylloflavan compound, −28.90 kcal/mol for Milk thistle, −31.83 kcal/mol for Ilexin B and −34.97 kcal/mol for Isosilybin B, as shown in Table 5.

MD Simulation
Docking analysis was used to determine the best position for the ligand to bind strongly to the receptor. The interaction patterns of the top complex with the target spike-receptor binding protein were determined using the MD simulation. After the MD simulation, the RMSD, RMSF, SASA and SSE distribution analyses were performed.

Root Mean Square Deviation (RMSD)
The RMSD was calculated during 100 ns simulations for the top complex. The RMSD was used as a function of time for the best ligand (Phylloflavan) with the spike-receptor binding domain protein. The RMSD results revealed that the MD simulation was equilibrated between 1 Å and 1.8 Å (Figure 4). Although it showed stability up to 50 ns, after that it showed minor deviations up to 90 ns; following that, it showed stability again. The RMSD results signified that the spike-receptor binding domain did not undergo large conformational changes. The RMSD plot of ligands in the right y-axis direction suggested that all these ligands were stable during simulation, with respect to the protein binding pocket.

Energy Calculations
The binding free energy (deltaG bind) of inhibitory drugs with a strong potential to inhibit the activity of the SARS-CoV-2 spike protein was calculated using MMGBSA. To perform the MMGBSA energy estimates, docking complexes with high energy function scores were obtained. By minimizing protein-ligand complexes, salvation energy and surface area energy, total free binding energy was determined. The results revealed that the compound name with ID had a high level of stability, with a total binding free energy of −30.35 kcal/mol for the Phylloflavan compound, −28.90 kcal/mol for Milk thistle, −31.83 kcal/mol for Ilexin B and −34.97 kcal/mol for Isosilybin B, as shown in Table 5.

MD Simulation
Docking analysis was used to determine the best position for the ligand to bind strongly to the receptor. The interaction patterns of the top complex with the target spikereceptor binding protein were determined using the MD simulation. After the MD simulation, the RMSD, RMSF, SASA and SSE distribution analyses were performed.

Root Mean Square Deviation (RMSD)
The RMSD was calculated during 100 ns simulations for the top complex. The RMSD was used as a function of time for the best ligand (Phylloflavan) with the spike-receptor binding domain protein. The RMSD results revealed that the MD simulation was equilibrated between 1 Å and 1.8 Å (Figure 4). Although it showed stability up to 50 ns, after that it showed minor deviations up to 90 ns; following that, it showed stability again. The RMSD results signified that the spike-receptor binding domain did not undergo large conformational changes. The RMSD plot of ligands in the right y-axis direction suggested that all these ligands were stable during simulation, with respect to the protein binding pocket.

Root Mean Square Fluctuations
The RMSFs were calculated to examine residential stability and flexibility over 100 ns. Trajectory conformation variability can be modified to compute the RMSFs for individual atoms. Phylloflavan had an average RMSF of 1.2 for overall positions. At the C-terminal and N-terminal, there were significant differences in all trajectories. The complexes' total residual fluctuations are depicted in Figure 5.
Molecules 2022, 27, x FOR PEER REVIEW

Root Mean Square Fluctuations
The RMSFs were calculated to examine residential stability and flexibility ov ns. Trajectory conformation variability can be modified to compute the RMSFs fo vidual atoms. Phylloflavan had an average RMSF of 1.2 for overall positions. At terminal and N-terminal, there were significant differences in all trajectories. The plexes' total residual fluctuations are depicted in Figure 5. SASA stands for solvent accessible surface area, and it is a new approach to proteins stable and folded. Figure 6 shows the computed SASA values for wild typ mutants. The average SASA values for Phylloflavan were determined to be 40, 80, 1 160, respectively, indicating that the available area of all the systems did not vary over the simulation procedure.

Water Bridges, Ionic Interactions and Hydrogen Bonding Graphs
Understanding the binding mechanism of both complexes in the binding po the spike-receptor binding domain necessitates atomic-level knowledge. Hydrogen ing, ionic interaction, hydrophobic contacts and salt bridges are all essential interm lar interacting forces in binding mode analysis. Over 100 ns simulation studies, th termolecular forces of interactions were predicted (Figure 7).

Solvent Accessible Surface Area (SASA)
SASA stands for solvent accessible surface area, and it is a new approach to keep proteins stable and folded. Figure 6 shows the computed SASA values for wild types and mutants. The average SASA values for Phylloflavan were determined to be 40, 80, 120 and 160, respectively, indicating that the available area of all the systems did not vary much over the simulation procedure.

Root Mean Square Fluctuations
The RMSFs were calculated to examine residential stability and flexibility over 100 ns. Trajectory conformation variability can be modified to compute the RMSFs for individual atoms. Phylloflavan had an average RMSF of 1.2 for overall positions. At the Cterminal and N-terminal, there were significant differences in all trajectories. The complexes' total residual fluctuations are depicted in Figure 5. SASA stands for solvent accessible surface area, and it is a new approach to keep proteins stable and folded. Figure 6 shows the computed SASA values for wild types and mutants. The average SASA values for Phylloflavan were determined to be 40, 80, 120 and 160, respectively, indicating that the available area of all the systems did not vary much over the simulation procedure. Figure 6. Graphs of solvent accessible surface area (SASA) revealed the protein's stability and folding.

Water Bridges, Ionic Interactions and Hydrogen Bonding Graphs
Understanding the binding mechanism of both complexes in the binding pocket of the spike-receptor binding domain necessitates atomic-level knowledge. Hydrogen bonding, ionic interaction, hydrophobic contacts and salt bridges are all essential intermolecular interacting forces in binding mode analysis. Over 100 ns simulation studies, these intermolecular forces of interactions were predicted (Figure 7).

Water Bridges, Ionic Interactions and Hydrogen Bonding Graphs
Understanding the binding mechanism of both complexes in the binding pocket of the spike-receptor binding domain necessitates atomic-level knowledge. Hydrogen bonding, ionic interaction, hydrophobic contacts and salt bridges are all essential intermolecular interacting forces in binding mode analysis. Over 100 ns simulation studies, these intermolecular forces of interactions were predicted (Figure 7).  Figure 7. After MD simulation, strong water bridges, ionic contacts and generated hydrogen bonding graphs can be seen in docked Phylloflavan in complex with spike glycoprotein.

Discussion
The best method to grasp the host-specific, infectious and pathogenic nature of COVID-19 is to fully comprehend the process for identifying viral receptors, which is critical for the creation of remedial treatments, medicines and antiviral cures [30]. There is mounting evidence that ACE2 gene polymorphism affects the association between ACE2 and the SARS-CoV-2 spike protein, influencing the viral entrance into the host organism and contributing to COVID-19 lungs and systemic harm [31]. Until now, none of the medications or treatments were proven to be effective against SARS-CoV-2, and the use and development of new medications is expensive and time-consuming due to the numerous testing stages required. [32][33][34]. Worse, the WHO predicts that COVID-19 will become endemic, signaling to the scientific community that treatment development for it has grown more critical. In this sense, docking methods are technically developed as a highly specific inhibitor for viral proteins, as well as an antiviral medication research and development process [35]. The worldwide risk posed by COVID-19 has prompted a rapid search for a therapeutic agent and the use of biological technologies, such as docking studies, as well as bioinformatics tools, to quickly assess the efficiency of other available medications against severe acute respiratory 2 infections [36].
A spike protein is a type I glycoprotein that protrudes from the surface of the virus and is the first component to make contact with the host cell. The establishment of COVID-19 treatment techniques is critical since the spike protein is the virus' key component for binding receptors to the host surface. An S-protein in the SARS-CoV-2 envelope is responsible for the interaction between host cells and the ACE2 receiver via the RBD unit [37]. This is the first and most important step in a SARS-CoV-2 infection and the interaction with the host cells. In the fight against SARS-CoV-2 inflammation, an antiviral medication that targets S-RBD-ACE2, which can transfer viruses into a host cell, gives a quick result [38,39].
Scientists have lately used computational capabilities to estimate the likelihood of joining numerous molecules before estimating and developing them in the lab. Molecular docking is used to find the binding pattern of small molecules against their target. As a result, molecular docking appears to be a valuable strategy for creating and screening novel chemicals to combat devastating diseases. The molecular interactions of 15,000 phytochemicals with different flavonoid subgroups and S-RBD-ACE2 were investigated. Structure optimization of ligand molecules and energy reduction in S-RBD-ACE2 were performed prior to docking analysis.

Discussion
The best method to grasp the host-specific, infectious and pathogenic nature of COVID-19 is to fully comprehend the process for identifying viral receptors, which is critical for the creation of remedial treatments, medicines and antiviral cures [30]. There is mounting evidence that ACE2 gene polymorphism affects the association between ACE2 and the SARS-CoV-2 spike protein, influencing the viral entrance into the host organism and contributing to COVID-19 lungs and systemic harm [31]. Until now, none of the medications or treatments were proven to be effective against SARS-CoV-2, and the use and development of new medications is expensive and time-consuming due to the numerous testing stages required [32][33][34]. Worse, the WHO predicts that COVID-19 will become endemic, signaling to the scientific community that treatment development for it has grown more critical. In this sense, docking methods are technically developed as a highly specific inhibitor for viral proteins, as well as an antiviral medication research and development process [35]. The worldwide risk posed by COVID-19 has prompted a rapid search for a therapeutic agent and the use of biological technologies, such as docking studies, as well as bioinformatics tools, to quickly assess the efficiency of other available medications against severe acute respiratory 2 infections [36].
A spike protein is a type I glycoprotein that protrudes from the surface of the virus and is the first component to make contact with the host cell. The establishment of COVID-19 treatment techniques is critical since the spike protein is the virus' key component for binding receptors to the host surface. An S-protein in the SARS-CoV-2 envelope is responsible for the interaction between host cells and the ACE2 receiver via the RBD unit [37]. This is the first and most important step in a SARS-CoV-2 infection and the interaction with the host cells. In the fight against SARS-CoV-2 inflammation, an antiviral medication that targets S-RBD-ACE2, which can transfer viruses into a host cell, gives a quick result [38,39].
Scientists have lately used computational capabilities to estimate the likelihood of joining numerous molecules before estimating and developing them in the lab. Molecular docking is used to find the binding pattern of small molecules against their target. As a result, molecular docking appears to be a valuable strategy for creating and screening novel chemicals to combat devastating diseases. The molecular interactions of 15,000 phytochemicals with different flavonoid subgroups and S-RBD-ACE2 were investigated. Structure optimization of ligand molecules and energy reduction in S-RBD-ACE2 were performed prior to docking analysis.

Data Collection and Ligand Database
For molecular docking, the 2D conformation of 15,000 phytochemicals were retrieved from six different databases, i.e., from MAPS database [40], PubChem [41], Zinc database [42], MPD3 [43], ChEMBL [44] and NPACT [45]. Each molecule of ligand was saved into the MOE database in mol form after minimization of energy. All of these were conducted through MOE. Drug designing, protein structure analysis, data processing and docking; all these steps were conducted through MOE.

Receptor Preparation and Analysis of Target Active Binding Sites
MOE performed molecular docking of the terpenoids dataset [46]. The l structure of the spike protein coupled to the ACE2 protein structure was retrieved from the PDB for docking purposes (PDB ID 6LZG) [47]. The MOE software's site finder function located the target protein's active sites. The ligand coordinates defined the active site in the original target protein site. The binding pocket containing the catalytic triad was selected.

Molecular Docking
Molecular docking is a computational procedure that is used to check the interaction between molecule and protein at an atomic level [48]. Various types of computational approaches are utilized to predict the affinity between the target and the active site [49]. Rigid ligand docking, semi-flexible ligand docking, and flexi docking are the three types. A particular scoring location confirms the ligand's approval in sub-atomic docking. There are specific scoring positions: in molecular docking, either a systematic or a ligand arch is used to evaluate a conformational search [50]. In molecular docking, binding energetics were tested using a force-field based scoring method, an empirical scoring function and a knowledge-based scoring function [51].

Analysis of Ligand Receptor Interaction
MOE's LigX tool was used to analyze the receptor-ligand interaction on 2D plots, in order to obtain a high-quality image of the top-docked complexes' receptor-ligand interaction. It showed a 2D graph of electrostatic contacts, hydrogen bonding, hydrophobic interactions and Van der Waals forces that were taken into account when determining the affinity of the drug-like molecule within the actively docked pockets. MOE is a program that creates three-dimensional pictures of protein-inhibitor complexes [52].

Physiochemical Property Profile and Toxicity Prediction
All the selected molecules were analyzed through the Molinspiration server for Lipinski's Rule of Five. This rule explains different drugs' properties, such as absorption, metabolism and drug secretion in the human body. Although the different drugs are evaluated on the basis of this rule, it also includes different values, such as MW, hydrogen HBA, HBD and log p-value. These are the standard of these values (A log p < 5, fewer than 10 H-bond acceptors, fewer than 5 H-bond donors and a molecular weight of fewer than 500 Daltons) [53]. Swiss-ADME software was used to predict the pharmacokinetic properties (absorption, metabolism, distribution, excretion and toxicity) [54].

MM-GBSA Binding Free Energy Calculations
Docking complexes were further validated by calculating the binding free energy (Prime/MM-GBSA) using the Schrodinger Suit Release 2020 [55]. The best poses of the inhibitory compounds related to the SARS-CoV-2 spike protein were chosen to obtain binding-free energies. The local optimization feature was used to minimize the docked complexes in Prime. Prime MM-GBSA was used to calculate the binding free energies based on the combination of the OPL-SAA force field, EMM (Molecular Mechanics Energies), the GSGB (solvation model for the polar solvation), the nonpolar solvation term (GNP) made up of solvent accessible surface area (SASA) and Van der Waals interactions. The binding free energy calculations were conducted based on the following equation: ∆Gbind = GComplex − Gprotein + Gligand.

Molecular Dynamic Simulation
Stability and effectiveness of the docked complexes were elucidated by performing the 100 ns molecular dynamic simulation [56]. The simulation was carried out using the Desmond package, which included a solvent system and a force-field called Optimized Potentials for Liquid Simulations 3 (OPLS3) [57]. Water molecules were used to solvate the molecular system. The system was electrically neutralized using Na/Cl ions. Before the MD simulation, the system used a heating process to reduce its energy use. Complexes were also reduced utilizing a steepest descent steps-based approach as part of the reduction method. Furthermore, after 1,000 sharpest descent steps, the system was brought into balance. Finally, the simulation was performed for 100 ns time at 300 K temperature and at 1 atm pressure using the NPT-ensembles.

Conclusions
The current study used an in silico approach encompassing many phases to recognize and identify active pharmacological targets for intimate drug creation against SARS-CoV-2. The goal of this research was to find druggable targets. Disease searching, disease-related genes, docking interaction and evaluation of docked complexes for medication potency were some of the bioinformatics and computational techniques utilized to follow the workflow. To select the compounds having the best residue interaction with the target protein, molecular docking was used. After examining the top 100 docked compounds, six molecules were chosen as the best molecules based on their docking score and drug evaluation. Phylloflavan, Milk thistle, Ilexin B and Isosilybin B were found as promising phytochemicals in the current investigation, with a good binding capability to SARS-CoV-2 S-RBD, and all have drug-like qualities. The findings of this study can be used to create and develop new drugs that have improved inhibitory activity on the spike glycoprotein. Although, for findings validation, an in vitro and in vivo experimental study is recommended.

Data Availability Statement:
The data generated and analyzed during this study are included in this article.