In Silico Identification and Analysis of Potentially Bioactive Antiviral Phytochemicals against SARS-CoV-2: A Molecular Docking and Dynamics Simulation Approach

SARS-CoV-2, a deadly coronavirus sparked COVID-19 pandemic around the globe. With an increased mutation rate, this infectious agent is highly transmissible inducing an escalated rate of infections and death everywhere. Hence, the discovery of a viable antiviral therapy option is urgent. Computational approaches have offered a revolutionary framework to identify novel antimicrobial treatment regimens and allow a quicker, cost-effective, and productive conversion into the health center by evaluating preliminary and safety investigations. The primary purpose of this research was to find plausible plant-derived antiviral small molecules to halt the viral entrance into individuals by clogging the adherence of Spike protein with human ACE2 receptor and to suppress their genome replication by obstructing the activity of Nsp3 (Nonstructural protein 3) and 3CLpro (main protease). An in-house library of 1163 phytochemicals were selected from the NPASS and PubChem databases for downstream analysis. Preliminary analysis with SwissADME and pkCSM revealed 149 finest small molecules from the large dataset. Virtual screening using the molecular docking scoring and the MM-GBSA data analysis revealed that three candidate ligands CHEMBL503 (Lovastatin), CHEMBL490355 (Sulfuretin), and CHEMBL4216332 (Grayanoside A) successfully formed docked complex within the active site of human ACE2 receptor, Nsp3, and 3CLpro, respectively. Dual method molecular dynamics (MD) simulation and post-MD MM-GBSA further confirmed efficient binding and stable interaction between the ligands and target proteins. Furthermore, biological activity spectra and molecular target analysis revealed that all three preselected phytochemicals were biologically active and safe for human use. Throughout the adopted methodology, all three therapeutic candidates significantly outperformed the control drugs (Molnupiravir and Paxlovid). Finally, our research implies that these SARS-CoV-2 protein antagonists might be viable therapeutic options. At the same time, enough wet lab evaluations would be needed to ensure the therapeutic potency of the recommended drug candidates for SARS-CoV-2.


Introduction
The World Health Organization (WHO) reported that the prevalence of SARS-CoV-2 is spreading at an alarming rate, posing severe health problems. The most recent outburst of second wave of SARS-CoV-2 has turned into a worldwide catastrophe. Following the coronavirus (CoV) epidemic in China in December 2019, WHO classified SARS-CoV-2, as the newest candidate of the Coronaviridae family within Nidovirales order [1]. As of 3rd May 2023, WHO has received a report from around 765,222,932 diagnosed COVID-19 infections worldwide, with 6,921,614 fatalities [2,3]. According to available information, the virus can be transmitted often by close, indirect, or direct exposure to infectious persons, as well as contaminated secretions such as nasal droplets and saliva, and respiratory secretions released when an infected individual sneezes, coughs, or speaks [3]. It has been linked to a wide range of signs and symptoms, consisting of minor to severe illness, which varies from patient to patient. Complications might appear anywhere from two to fourteen days after the virus has been infected. Fever, fatigue, chronic cough, sore throat, difficulty breathing, impairment of taste/odor, nausea, sputum production, headache, expectoration, diarrhea, anorexia, and some other symptoms might occur at separate phases of the disease [1,4].
SARS-CoV-2 is a membrane-encased positive-sense single-stranded RNA ((+) ssRNA) virus having a diameter ranging from 60 to 140 nanometers [4,5]. The envelope is surrounded by spike-shaped glycoprotein protrusions that resemble crowns under the electron microscope [6]. The spike (S), nucleocapsid (N), envelope (E), and membrane (M) proteins are among the four crucial targets encoded by the SARS-CoV-2 genome. Main protease (3CLpro), RNAdependent RNA polymerase (RdRp), and papain-like protease (PLpro) are some of the nonstructural proteins synthesized by the viral DNA [7]. Nonstructural protein 3 (Nsp3) proteins containing macrodomains are pervasive and evolutionarily conserved and responsible for the transcription process [8]. Previous study has established that human angiotensin converting enzyme 2 (ACE2) receptor has a greater affinity for the RBD region of the spike protein [9]. The attachment of favorable ligands to the active pockets of human ACE2 receptor might alter the protein's structure. As a result, the viral ACE2 entrance region might be a feasible object for therapeutic advancement. Since the main protease of SARS-CoV-2 is vital for its growth and the consequent expression of the replicase polyproteins, it has turned into an obvious target for anti-COVID-19 therapeutic design [10]. As a result, focusing on these proteins might help with long-term COVID-19 infection management and eradication.
The viral disease is spreading at a surprising pace worldwide, and researchers are racing to develop effective drugs to use as therapeutic agents. The most promising choices appear to be natural compounds with substantial bioavailability and minimal cytotoxicity [1]. Clinically approved antiviral drugs are effective; however, some people become resistant to drugs. In contrast, it has been claimed that phytochemicals have more acceptable side effects and can be a satisfactory substitute for synthetic antiviral compounds for the suppression of viral life-cycle and penetration [10].
Humans have always relied on natural compounds, especially phytochemicals, to treat health problems since the dawn of time. Recently, Shawan et al. presented luteolin and abyssinone II as possible phytochemicals against SARS-CoV-2 [1]. Besides, Manojkumar et al. reported ervoside had anticoronavirus properties [11]. Similarly, Emran et al. identified phytochemicals medicagol, faradiol, and flavanthrin as the potential barrier of SARS-CoV-2 [12]. Computer-assisted drug development (CAD) entails the usage of computerized techniques to discover, design, and evaluate therapeutics and associated pharmacologically active substances [13]. CAD techniques have improved compound screening significantly over time, aimed at targeting structure prediction and model development, active site determination, comprehending the protein-ligand complex, testing a huge dataset of substances by estimating their pharmacokinetics characteristics, and analyzing the dynamics of proteins binding with ligands within biological settings [14]. Existing medicines like Molnupiravir and Paxlovid have been authorized by the FDA for utilization in emergencies; the treatment may be used either alone or combined with others [15]. For COVID-19 patients, the antiviral drug Molnupiravir has been recommended as a therapeutic for SARS-CoV-2 because it increases the likelihood of viral RNA alterations while also inhibiting viral replication [16]. Through inhibition of proteasome breakdown of viral proteins, Paxlovid inhibits protein production (RNA-dependent RNA polymerase, helicase, exoribonuclease, RNA-binding protein endoribonuclease). Consequently, the viral transcription and replication are halted [7].
The main focus of this in silico work was to utilize computational tools, i.e., molecular docking and MD simulation to examine the effective binding interactivity and affinities of repurposed antiviral phytochemicals with the human ACE2 receptor, Nsp3 macrodomain, and the main protease of the SARS-CoV-2 virus and identify the finest ligand hit [17]. Among all other crucial characteristics, absorption, distribution, metabolism, toxicity, and excretion (ADMET) were evaluated, and the best of them were selected. Finally, the most effective phytochemicals with higher binding energy to the target receptor and stronger stabilizing capacity were confirmed by employing molecular dynamics simulation.

Materials and Methods
Virtual screening of natural bioactive molecules has become the standard method in the present therapeutic development workflow [18]. In this study, a wide range of repurposed phytochemicals were used from the NPASS (http://bidd .group/NPASS/) and PubChem (https://pubchem.ncbi.nlm .nih.gov/) servers as prospective ligands for SARS-CoV-2. The recently approved COVID-19 antiviral drugs Molnupiravir and Paxlovid were used as control drugs [19]. The workflow of our work was provided in Figure 1.

Characterization of Drug-Likeness Properties.
A druglike molecule can be considered a drug candidate by assessing its drug-like properties. The canonical SMILE sequence of the 1163 small molecules was fetched from the PubChem drug web server. The free accessible SwissADME was employed to compute the major physicochemical descriptors, pharmacokinetic properties, drug-like parameters, and associated factors [20]. To analyze the results, this application employs five principles of Lipinski's rule [21], Ghose's rule [22], Veber's rule [23], Egan rule [24], Muegge's rule [24], the number of rotatable bonds, and TPSA. Swiss-PdbViewer was subsequently used to minimize the energy of the selected proteins [28]. Next, the energyminimized structures were loaded into AutoDock-MGLTools 1.5.6 to incorporate polar hydrogen and convert the PDB to PDBQT format.

Active Site Detection and Grid Box Preparation.
Finding a ligand-binding region on a protein is the basic strategy for the molecular docking technique [29]. The possibility of protein-ligand attachment relies on numerous factors such as hydrogen bonds, hydrophobic or hydrophilic interactions, electrostatic and salt bridges. CASTp 3.0 website (http://sts.bioe.uic.edu/castp/) was employed to detect the active region of target proteins [30]. It applies an alpha shape detection approach to determine topographic properties and estimate protein area and volume for identifying ligand-binding cavities.

Binding Affinity Prediction by AutoDock vina.
Virtual screening via docking studies is extensively used in computer-led pharmaceutical research to uncover promising drug-like substances. Initially, AutoDock vina was exploited to conduct rigid molecular docking among the proteins and selected compounds (ligands and control drugs), including a search area of 27,000 m 3 and exhaustiveness 10, and ligands being flexible while receptors remained rigid [18]. AutoDock vina calculates the binding energy and fixes the binding poses using the Lamarckian genetic algorithm. Here, in this study, 149 small molecules were docked with three target proteins (coordinates of geometry-optimized ligands of the best hits provided in Supplementary

BioMed Research International
Previously screened ligands having higher affinity for target proteins than the reference drugs were explored in this step.

Preparation of Ligand
Structures. The LigPrep module yields top hits of 3D configurations for small molecules, beginning from 1-dimensional/2-dimensional/3-dimensional structures in Maestro, Mol2, SMILES, or SD format [31]. By introducing hydrogens, ionizing at pH (7 ± 2:0), and subtracting salts, the LigPrep tool builds molecules and constructs 3D structures of them. Following that energy minimized and geometrically refined ligands were prepared by employing a built-in OPLS3e force field in Schrödinger Maestro 12.5 [32].

Preparation of Protein
Structures. The protein structures (main protease, Nsp3, human ACE2 receptor) were loaded straight into the protein preparation wizard [32]. Protein structures were preprocessed by setting up bond orders, adding hydrogens and cap termini, and filling the missing atoms by prime module. At pH 7.0, the PROPKA application was used to calculate the protonation phases. Following that, the water portion around the protein was eliminated above 3.0 Å, and restrained minimization was executed utilizing the OPLS3e force field.

Preparation of Receptor Grid
Box. The grid region directs small molecules to the binding center of the protein, making it an important part of molecular docking research. The grid model was created with the standard options of a Van der Waals radius scaling marker of 1.0 and a charge threshold score of 0.25 in the receptor grid generation package. The attached ligands UQZ, AR6, and XX5 within the protein structures main protease, Nsp3, and human ACE2 receptor, respectively, were used to define the region for the grid map.

Glide Docking and MM-GBSA Studies.
Glide is a combined molecular docking technology that can facilitate both ligand and receptor flexibility [33]. Glide XP was developed to retrieve the finest docking poses having the greatestscoring compounds. For drug molecules, a minimum scoring of 0.15 and a Van der Waals radius scaling marker of 0.80 was applied.
Docking score = a × VdW + b × Coul + Hbond + Metal Here, a and b are coefficient constants for VdW and Coul, respectively, VdW is the Van der Waals energy, Coul is the Coulomb energy, H-bond is the hydrogen bonding with the receptor, Metal is the binding with metal, Lipo is the constant term for lipophilic, BuryP is the buried polar group penalty, RotB is the rotatable bond penalty, and Site is the active site polar interaction [1].
The binding free energy among the receptor and the collection of small molecules was measured using the prime MM-GBSA module. The binding energy of the ligandprotein constructs was estimated utilizing the OPLS3e force field, and the docked conformations were minimized utilizing Prime's native optimization tool.
ΔGbind Binding Free Energy ð Þ = ΔEMM + ΔGsolv + ΔGSA: Here, ΔEMM represents lowered energy deviations among the totality of the energies of the protein and ligand and protein-ligand complex. ΔGsolv displays the divergence in the GBSA solvation energy of the complex structure and the aggregate of the salvation energies for the ligand and protein. ΔGSA describes the deviation in the energies for the surface area of a complex and the total surface area of the ligand and protein complex [34].
2.5. Molecular Dynamics Simulation by GROMACS. The molecular dynamics program simulates the movements of a protein molecule utilizing the interaction potential to compute interatomic energies and equations of motion that regulate the machinery's dynamics in the drug design study. It illustrates the stability and flexibility data of ligand binding to a flexible target protein. GROMACS (https://simlab .uams.edu/) service was exploited to simulate the proteinligand conformations, and the GROMOS96 43a1 force field was employed to produce the topological data of the complex constructs [35]. The PRODRG (http://davapc1.bioch .dundee.ac.uk/cgi-bin/prodrg) Server was employed to render small molecule topology and coordinate information [36]. The aqueous phase of macromolecules was produced sequentially using the SPC water model (simple pointcharge) and subsequently neutralized using 0.15 M NaCl solution [37]. A triclinic box was used to contain the bimolecular environment, and 5000 iterations of steepest descent strategies were used to minimize energy. The equilibrium of ion molecules around the macromolecule was accomplished at 310 K and 1.0 bar utilizing NPT (constant pressure) and NVT (constant volume) setups. After 100 nanoseconds of simulation, it provided trajectories of simulated structures, including the root-mean-square deviation (RMSD), the rootmean-square fluctuation (RMSF), the solvent-accessible surface area (SASA), hydrogen bonds (HBs), and the radius of gyration (Rg) [38].

Molecular Dynamics Simulation and Post MM-GBSA
Evaluation by Desmond. The molecular dynamics simulation provides evidence regarding the mobility and stability of the bound protein-ligand complex. On Desmond software, the MD simulation and post-MMGBSA analysis of the main-protease_ligand, Nsp3_ligand, and human ACE2 receptor_ligand complexes were performed [39]. These compounds were solvated on a cubic TIP3P water model using the system builder package. A minimal spacing of 10 was maintained between the protein and the solvated region. Subsequently, Na+ salts were supplied until the final system strength reached 0.15 M, which is the physiological salt concentration present in the human body. The integrated OPLS3e force field was used to optimize the final system's energy. To complete the MDS, we used the isothermal isobaric ensemble (NPT) at 1.013 bar and 310°K. The total

12
BioMed Research International duration of the simulation run was 100 nanoseconds (ns). It was paired with a recording duration of 100 picoseconds (ps), during which 1000 frames were incorporated into the trajectory. Next, we studied the trajectories in the simulation interaction diagram (SID) program, and the reported results comprised RMSD, RMSF, protein-ligand contact outline, and biophysical properties of ligands. After running the simulations, MM-GBSA was evaluated employing the thermal MM-GBSA.py program. During the assessment, a 0-1000 periodic frame was incorporated for the analysis [40].

Prediction of Molecular Target with SwissTargetPrediction
Server. The anticipation of a molecular target for a smallmolecule is vital for drug research and development. These studies are essential for assessing the potential for adverse reactions or cross-reactivity in Homo sapiens caused by the action of that bioactive small molecule. We employed Swis-sTargetPredcition (http://www.swisstargetprediction.ch/) to determine the human body receptors for small compounds that had previously been identified by molecular docking and shown stability via MD simulation investigations [41]. The canonical smiles of the small compounds were used in the server and analyzed.

Prediction of Biological Activity by PASS-Way2Drug Tool.
The PASS-Way2Drug webserver (http://www.pharmaexpert .ru/passonline/) was employed to the prediction of biological activity scales for Lovastatin, Sulfuretin, and Grayanoside A [42]. For PASS recommendations to be reliable, the Pa (likelihood to be effective) threshold should be set at 70% or above. This is because surpassing the Pa>70% threshold yields very reliable predictions [42]. Calculated ligand activity was based on Pi and Pa scores.          BioMed Research International shortlisted for the following evaluation (Supplementary table1). Table 1 represented the drug-like properties of the best-hit phytochemicals and control drugs.

Analysis of ADMET Properties.
A total of 149 drug-like substances were qualified after this analysis. Moreover, from estimating distribution levels, all the compounds are impermeable to the blood-brain barrier. Metabolic inability could cause lower bioavailability and excretion, high toxicity, and drug-drug interactions. These 149 small substances function as isoforms of the CYP 2D6 and 3A4 enzymes. Diverse computational algorithms are used to evaluate toxicity: hERG inhibitors, AMES toxicity, maximum tolerated dosage Hepatotoxicity. Ligands with a negative value in these models were chosen for the following step. ADMET properties of the best-hit phytochemicals and control drugs were presented in Table 2. Finally, we filtered out 149 drug-like substances from this analysis (Supplementary table2).  Table 3c). Lovastatin's binding energy for the main protease was -7.2 kcal/mol, which was considerably higher than that of the control ligands, Molnupiravir (-6.4 kcal/mol), and Paxlovid (-6.6 kcal/mol) ( Table 3). Lovastatin produced a robust hydrogen interaction with the amino acids ARG131 (2.30102 Å) of the main protease, whereas Molnupiravir and Paxlovid formed three and six hydrogen bonds with the target protein, respectively, with THR26, HIS41, ASN119, ASN142, GLY143 (1.98365 Å), and CYS145 residues (AutoDock vina). Sulfuretin had binding energy of -8.8 kcal/mol for Nsp3 compared to the control ligands molnupiravir and Paxlovid, which had binding energies of -7.7 and -7.5 kcal/mol, respectively (Table 4). Sulfuretin formed seven strong hydrogen bonds with the Nsp3 protein (VAL49, LEU126, SER128, ALA129, GLY130, PHE156, and ALA38), whereas Molnupiravir and Paxlovid created  Table 5). Molnupiravir and Paxlovid formed five (ASP206, HIS378, ASN394, ARG514, and LYS562 (2.198 Å)) and six (ASP206, ALA348, TRP349 (1.978 Å), ASP350, HIS378, and ARG514) hydrogen bonds with the target protein, human ACE2 receptor, respectively. Grayanoside A formed three strong hydrogen bonds (ARG273, ARG273, and GLU406) and six hydrophobic bonds (VAL209, LYS562, TRP566, LEU95, LYS562, and ALA99) with the human ACE2 receptor.   (Figures 2-4), and it showed the comparative representation of protein-ligand complexes of the best hit ligands and control drugs. Here, Tables 6-8 summarized the Glide score and MM-GBSA scores between three target proteins and the best hit phytochemicals and control drugs. Analysis of glide XP score and MMGBSA values, it was evident that Lovastatin is better candidate than other potential ligands. It formed three hydrogen bonds (HIS163,  5-7). Lovastatin, Sulfuretin, and Grayanoside A were found inside the binding cavity with the cocrystallized compound ( Figure 8). As a result, they were proved to be the best candidate for main protease and Nsp3 of SARS-CoV-2 and human ACE2 proteins, respectively.

Analysis of Molecular Dynamics Simulation.
To circumvent the fundamental drawback of molecular docking, we ran a computational MD simulation, which incorporated the dynamic character of the protein following inhibitor binding. This experiment produced statistical figures for the RMSD, RMSF, hydrogen bonds, SASA, and Rg values of protein-ligand complexes. The average RMSD of main protease_Lovastatin, main protease_Molnupiravir, and main protease_Paxlovid complexes for the main protease was 0.312696293 nm, 0.291836715 nm, and 0.326214306 nm, respectively, indicating that the chosen drug candidate Lovastatin exhibited an identical result compared to Molnupiravir and Paxlovid ( Figure 9). As per Figure 9, the main protease_Lovastatin and main protease_Molnupiravir complexes were stable with a fixed RMSD value less than 0.35 from 30 to 80 ns, but the main protease_Paxlovid complex had an increased RMSD value more than 0.35 after 75 ns. Similarly, the predicted average RMSD values of the ligands  The average RMSD value of the Nsp3_Sulfuretin and Nsp3_Paxlovid complexes for SARS-CoV-2 Nsp3 protein was 0.297815 nm and 0.284552759 nm, respectively, though the Nsp3_Molnupiravir complex displayed an increased variation of RMSD value exceeding 0.35 nm after 45 ns ( Figure 10). During the 100 ns simulation timeline with Nsp3 protein, control drugs molnupiravir and Paxlovid The RMSD value of the human ACE2 receptor_Grayanoside A, human ACE2 receptor_Molnupiravir, and human ACE2 receptor_Paxlovid complexes for human ACE2 protein stayed under 0.35 nm, and stable throughout the 100 ns run ( Figure 11). Likewise, the ligands Grayanoside A, Molnupiravir, and Paxlovid had average RMSD values of 0.601344 nm, 0.933326 nm, and 0.43800 nm, respectively. The RMSF of backbone heteroatoms per residue of the human ACE2 recep-tor_Grayanoside A complex stayed within 0.4 nm, with higher RMSF oscillation inside loops from 137 to 139 and 338 to 340 residues. On the other hand, peaks inside loop regions beyond 0.4 nm were evident from 137 to 140 and 331 to 345 residues for human ACE2 receptor_Molnupiravir and human ACE2 receptor_Paxlovid complexes, respectively. The average Rg values of human ACE2 receptor_Grayanoside A, human ACE2 receptor_Molnipiravir, and human ACE2 receptor_ Paxlovid complexes were 2.329435, 2.342172667, and 2.335405325 nm, respectively. The average hydrogen bond interactions for the complexes human ACE2 receptor_Grayanoside A and human ACE2 receptor_Molnupiravir were 499.0 and 492.0, respectively, whereas the complex human ACE2 receptor_Paxlovid had a higher value of 498.0. Figure 11 shows that the SASA values of the human ACE2 receptor_Grayanoside A, human ACE2 receptor_Molnupiravir, and human ACE2 receptor_Paxlovid complexes were stable throughout the simulation, suggesting that the protein's surface area remained unaltered.  (Figure 12). Sulfuretin interacted with 7NT3 creating bonds with THR26 (hydrogen bonds and water bridges), GLY143 (hydrogen bonds and water bridges), SER144 (hydrogen bonds and water bridges), CYS145 (hydrogen bonds and water bridges), and GLU166 (hydrogen bonds, water bridges, and ionic bonds) amino acids for 30%, 20%, 30%, 40%, 20%, and 100% of 100 ns timeframe. Molnupiravir interacted with HIS41 (hydrophobic), GLU166 (water bridges), VAL186 (hydrogen bonds and water bridges), and GLN189 (hydrogen bonds and water bridges) for 80%, 100%, 70%, and 90% of 100 ns timescale. Paxlovid had bonds with HIS41 (hydrophobic, hydro-gen bonds, and water bridges), GLY143 (hydrogen bonds and water bridges), SER144 (hydrogen bonds and water bridges), and GLU166 (hydrogen bonds, water bridges, and ionic bonds) for 50%, 90%, 40%, and 300% of the simulation period ( Figure 13).
3.8. Analysis of Activity Spectra of the Phytochemicals. Using the identified compounds, prediction of activity spectra for substances (PASS) was carried out and is shown in Supplementary Tables 6a, 6b, 6c. In our study, we used PASS to build a predictive model, and we kept the Pa (likelihood of activity) that was higher than 70%; since an absolutely durable forecast may be made using the Pa > 70% (0.7) criteria. Lovastatin had 18 biological activities, Sulfuretin showed 27 activities, and Grayanoside A possessed 30 biological features. Lovastatin, Sulfuretin, and Grayanoside A present Pa values greater than 0.70 across the board to be considered for use as an active biological agent in the treatment of SARS-CoV-2.

Discussion
In recent years, pandemics and epidemics caused by viruses have become one of the most prevalent reasons for infections and mortality worldwide. SARS-CoV-2, the updated variant of coronaviruses, has led to a catastrophe, with 665,518,891 and 6,714,212 confirmed cases and fatalities, respectively (11th January,2023; https://covid19.who.int/). Surprisingly, currently, a limited amount of effective anti-SARS-CoV-2 therapeutics are available, and most of them are under investigation [43].
Following the outbreak of SARS-CoV-2, Mpro, also regarded as 3CLpro (main protease), became a viable therapeutic focus due to its involvement in the development of replication-translation mechanisms. Furthermore, given the accessibility of high-resolution protein structures, these proteins have a highly conserved sequence and no homology with human proteins [44]. Nsp3 is a multidomain protein with a Glu-rich acidic domain at the N-terminus, an X domain, a SUD domain, a papain-like protease domain, and a transmembrane domain. Nsp3 is responsible for viral multiplication and pathogenesis in humans and facilitates immune evasion via its hydrolyzing capability [43].  LYS_28  ALA_38  ALA_39  ASN_40  LYS_44  GLY_46  GLY_47  GLY_48  VAL_49  ALA_50  ALA_56  THR_57  ASN_58  ALA_60  GLN_62  VAL_63  GLU_64  ASP_66  ASP_67  ALA_70  SER_84  HIS_86  ASN_87  ASN_99  VAL_100  ASN_101  LYS_102  GLY_130  GLY_133  ILE_131  PHE_132  PHE_156  ALA_134  ASP_157  ASN_159  LEU_160

24
BioMed Research International ACE2 receptor on the cellular wall permits the virus to enter cells, which is required for infection to begin [45]. To inhibit these viral proteins, we utilized phytochemicals with druglike properties. This research was divided into three sections, namely, virtual screening (VS) of the physicochemical and pharmacokinetic features of drug-like compounds, virtual screening by molecular docking of proteins and ligands, and simulation of the best hit complexes. In the initial stage, we studied the drug-like characteristics of the ligands utilizing the five principles of Lipinski. We stuck to the established guidelines, i.e., hydrogen bond donors ≤ 5, hydrogen bond acceptors ≤ 10, molecular mass < 500, and logP < 5. The molecular weight of a small molecule can influence its absorption, bile excretion ratio, blood-brain barrier passage, and engagements involving receptors [46]. Likewise, hydrogen donor and hydrogen acceptor groups are mostly responsible for the permeability and polarity of the drug-like molecules. Lipophilicity is an indicator that influences the metabolism and solubility of those molecules. A lower or higher score might impede this characteristic [47]. TPSA refers to the area belonging to polar atoms like nitrogen, oxygen, and their associated hydrogens [48]. Out of 1163 small molecules, 497 of them passed the criteria. We tested the pharma-cokinetic figure of the ligands before analyzing their binding affinity and orientation. The characteristics of a small molecule in terms of ADMET properties make it a viable candidate. Using the human intestinal absorption (HI) prediction score and the Caco-2 permeable theory, the probability that the small molecules would reach systemic circulation and exert their function was calculated [49]. Pglycoprotein serves as a drug carrier and eliminating compounds from different organs [50]. The main subfamily (2D6, 2C9, and 3A4) of cytochrome P450 monooxygenase (CYP) enzymes is crucial in the metabolism of the druglike molecules [51]. In the initial phases of pharmaceutical research, AMES mutagenicity is commonly used to determine the probability of genotoxicity and teratogenicity [52]. Cardiovascular poisoning might be caused by inhibiting the cardiac human ether-a-go-go-related (hERG) gene [53]. We also tested the maximum tolerated dose of chemical substances for the human body. Eventually, only 149 drug-like molecules passed the ADMET evaluation.
The importance of virtual screening employing molecular docking has grown significantly in the field of drug development over time. According to the study, 24 small molecules had a greater binding affinity against the main protease (7NT3) than the reference compounds:  TYR_127  LEU_144  GLU_145  ASN_149  ASP_269  ARG_273  PHE_274  THR_276  VAL_343  CYS_344  HIS_345  PRO_346  THR_347  ALA_348  LYS_363  ASP_367  ASP_368  THR_371  HIS_374  GLU_375  HIS_378  GLU_402  GLU_406  THR_445  THR_449  LEU_503  PHE_504  HIS_505  TYR_510  SER_511  PHE_512  ARG_514  ARG_518 SER_44  SER_47  ASP_206  ARG_273  HIS_345  PRO_346  THR_347  ALA_348  TRP_349  ASP_350  HIS_374  GLU_375  HIS_378  ASP_382  ASN_394  GLU_398  HIS_401  GLU_402  PHE_504  HIS_505  ASN_508  ASP_509  TYR_510  SER_511  ARG_514  TYR_515 H-Bonds Hydrophobic  LYS_363  THR_365  ASP_367  ASP_368  LEU_370  THR_371  HIS_374  GLU_375  HIS_378  GLU_402  GLU_406  SER_409  LYS_441  GLN_442  THR_445  SER_502  PHE_504  HIS_505  TYR_510  ARG_514  TYR_515 [55]. It was also found inside the binding pocket (Asp22, Ile23, Gly48, Val49, Gly130, or Phe156) [43]. A total of 20 phytochemicals had higher glide XP scores over 6 kcal/mol and 2 of them showed binding-free energy above -50 kcal/mol against human ACE2 receptor. Based on the glide and MMGBSA scores, as well as the number of hydrogen bonds, we selected Grayanoside A as the lead candidate against human ACE2 receptor (1R4L). Most of the residues of the active site (Tyr515, Arg514, His505, Phe504, Glu402, His378, Glu375, His374, Asp368, Cys361, His345, Cys344, and Glu145) were found attached to Grayanoside A [56]. Pai et al. found that isochlorogenic acid showed inhibition activity against human ACE2 receptor with a glide score MM-GBSA values of −8.799 and −44.248. MD simulations offer a plethora of energetic data on protein and ligand binding, as well as a wealth of structural figures on biological macromolecules. This type of knowledge is crucial for comprehending the structural and functional correlation of the receptor and the basis of proteinligand association, and also for steering the therapeutic research [51]. During the simulation trajectory, the RMSD of the protein Cα and RMSF of the amino acids, also the ligand-protein H-bonding association, the solventaccessible surface area (SASA), and the radius of gyration (Rg), were assessed to determine the steadiness of the complex structures [52]. The RMSD value is considered to indicate the flexibility and dynamic character of the protein. It showcased the movement of amino acids along with the MD simulation [38]. Thus, a relatively large RMSD value suggested more motion, whereas a relatively low RMSD value of protein showed less movement. The RMSD results suggested that the RMSD value of the main protease_Lovastatin backbone was identical to those of the reference complexes main protease_Molnupiravir and main protease_ Paxlovid. The ligands Lovastatin and Paxlovid remained steady, with two short peaks. As a result, the protein might remain stable during the simulation, after the attachment of the Lovastatin molecule. A detailed investigation of the RMSF demonstrated the specific fluctuation of amino acids in the catalytic and noncatalytic areas of the protein-ligand complexes. The RMSF value confirmed that the attachment of Lovastatin to the receptor might not increase flexibility.
The Rg values display the compactness of the protein with folding and unfolding nature by the thermodynamic concept. The interaction of the ligand Lovastatin did not modify the structure of the protein. Hydrogen bonds are another vital determinant of protein-ligand stability. Protein-ligand interaction is stronger with more hydrogen bonds. When compared to the reference complexes, the main protease_Lovastatin complex had a similar amount of hydrogen bonds, indicating a stable protein-ligand construct. The unfolding of the protein during the denaturation process exposes nonpolar hydrophobic interactions to the aqueous system. As a result, the protein's structure is disrupted. The SASA value computing determines the fluctuation in the accessibility of protein to solvent [57]. The SASA analysis revealed a similar tendency, with both main protease_Lovastatin and main protease_Molnupiravir complexes exhibiting significant similarities throughout the simulation.     The RMSD values suggested that Lovastatin, Sulfuretin, and Grayanoside A were firmly bound to the proteins. The RMSF plots implied that the main protease (3CLpro), Nsp3, and human ACE2 receptor were structurally stable while bound with respective ligands. The protein-ligand attachment maps continuously showed that the proposed ligands-maintained contact with active site residues. Lastly, the postsimulation MM-GBSA results stated that Lovastatin, Sulfuretin, and Grayanoside A had a higher free-binding affinity for their respective proteins. Previous structure-based computational research yielded similar findings, demanding wet-lab investigation. According to study lead by Gurung et al., bonducellpin D was found as a potential inhibitor for SARS-CoV-2 3CLpro protein [58]. In another study, Eissa et al. identified vidarabine as prospective antiviral agent for SARS-Cov-2 nonstructural protein-10 [59]. Ottavia Spiga et al. found simeprevir and lumacaftor the most potent inhibitors of Spike protein on the basis of their computational findings [60]. Kusumaningsih et al. found luteolin and naringenin as the probable drug candidates for main protease of SARS-CoV-2 [61]. Lovastatin, Sulfuretin, and Grayanoside A have been reported as antiviral agents [62][63][64]. Our structure-based strategy again showed antiviral activity of these small substances against SARS-CoV-2 critical protein targets. However, these compounds should be examined further in the pharmaceutical research facility to evaluate their potency, inhibitory power, and toxicity against their respective targets. While there is no denying the enormous success of drug repurposing, the in silico approach is not without its limitations. In a similar fashion, one disadvantage of molecular docking is the lack of proper scoring functions and algorithms, which may compromise molecular screening. Another challenge that researchers face is the difficulty in selecting the most effective target combinations due to the absence of quantifiable data for assessing network dynamics, as well as the inability to construct the molecular network of the disease [18,65]. Apart from those certain constraints due to data reliability, biasness and irregularities in the available current data, our study shows a comparison between established compounds and screened compounds using several bioinformatics tools and introduces in silico models that can swiftly present a summary of prospective therapeutic options economically and expediently for a microorganism such as SARS-CoV-2, which is constantly mutating and without any established therapy.

Conclusion
Repurposing drug-like phytochemicals is a secure means of building new therapeutics, with the main benefit of lowering the cost and duration of preclinical trials for novel candidates. The COVID-19 infectious disease induced by the SARS-CoV-2 outbreak has caused a worldwide medical catastrophe and finding a suitable cure for the virus continues to be a primary concern. The findings of our work indicated that using a structure-based strategy such as molecular docking and MD simulations, novel therapeutic candidates may be developed that selectively address the nonstructural protein 3, the main protease from SARS-CoV-2, and the human ACE2 protein. A preliminary screening of 1163 small phytochemicals combining drug-likeness and ADMET characteristics resulted in the identification of 149 of them. The degree of binding interaction and energy between the filtered compounds and the main protease, nonstructural protein 3, and human ACE2 receptor was estimated utilizing the docking procedure on the AutoDock vina and Schrodinger Suites. Compounds Lovastatin, Sulfuretin, and Grayanoside A outperformed Molnupiravir and Paxlovid in terms of binding score and hydrogen bond numbers against the main protease, Nsp3, and human ACE2 receptor, respectively. Eventually, 100 ns MD simulation studies of 3CLpro_ligand, Nsp3_ligand, Grayanoside A_ ligand complexes were completed to evaluate and improve our proposed design. This investigation is aimed at determining the promising inhibitors and devise protocols for continual improvement of COVID-19 medications. To summarize, all the repurposed compounds suggested here may provide a holistic understanding of structure-based drug development for SARS-CoV-2 given that they continue to remain potent in further drug development processes.

Data Availability
Availability of additional data (supplementary files) will be provided on request.

Conflicts of Interest
There are no conflicts of interest to declare.