Molecular Docking and Dynamics Simulation of Several Flavonoids Predict Cyanidin as an Effective Drug Candidate against SARS-CoV-2 Spike Protein

The in silico method has provided a versatile process of developing lead compounds from a large database in a short duration. Therefore, it is imperative to look for vaccinations and medications that can stop the havoc caused by SARS-CoV-2. The spike protein of SARS-CoV-2 is required for the viral entry into the host cells, hence inhibiting the virus from fusing and infecting the host. This study determined the binding interactions of 36 flavonoids along with two FDA-approved drugs against the spike protein receptor-binding domain of SARS-CoV-2 through molecular docking and molecular dynamics (MD) simulations. In addition, the molecular mechanics generalized Born surface area (MM/GBSA) approach was used to calculate the binding-free energy (BFE). Flavonoids were selected based on their in vitro assays on SARS-CoV and SARS-CoV-2. Our pharmacokinetics study revealed that cyanidin showed good drug-likeness, fulfilled Lipinski's rule of five, and conferred favorable toxicity parameters. Furthermore, MD simulations showed that cyanidin interacts with spike protein and alters the conformation and binding-free energy suited. Finally, an in vitro assay indicated that about 50% reduction in the binding of hACE2 with S1-RBD in the presence of cyanidin-containing red grapes crude extract was achieved at approximately 1.25 mg/mL. Hence, cyanidin may be a promising adjuvant medication for the SARS-CoV-2 spike protein based on in silico and in vitro research.


Background
As the COVID-19 pandemic is at the end of its third year, public health ofcials must assess where we are and how we may break the SARS-CoV-2 devastating grip in the world. Te rapid discovery of many safe and efective COVID-19 vaccines has been one of the pandemic's biggest scientifc achievements. However, vaccines alone will not be enough to stop the pandemic due to more transmissible new variants, and vaccines are designed to guard against severe sufering and death [1]. New variants due to mutation in SARS-CoV-2 are a subject of concern because the emerging mutants have a potential for enhanced infectivity, competitive ftness, and transmission [2]. Tere occurs an accumulation of mutations, which drives viral evolution and genome variability, and this enables the virus to escape host immunity and develop drug resistance [3]. Important mutations have appeared in the SARS-CoV-2 spike protein that interacts with the host immune system [4]. Te reported variant, omicron (B.1.1.529) [5], has 30 mutations in the spike protein along with K417N, which was earlier recognized to reduce the efectiveness of a cocktail of therapeutic monoclonal antibodies [6].
Te SARS-CoV-2 spike protein is a heavily glycosylated homotrimeric transmembrane protein; the S1 subunit contains an RBD that facilitates viral attachment to the host receptor angiotensin-converting enzyme-2 (ACE2), and the S2 subunit mediates viral entry by mediating host-viral membrane fusion [7]. Te biomechanical strength of the ACE2-spike protein manages viral adherence and access to host cells [8]. Hence, the spike protein of SARS-CoV-2, which plays a crucial role in viral attachment and fusion, could be a potential therapeutic target.
Several strategies have been used to develop antiviral drugs for SARS-CoV-2; to date, some drugs are efective against this virus, but only some neutralizing antibodies and remdesivir have been approved by the US FDA [9]. Te drug development against SARS-CoV-2 is mainly focused on interrupting the virus's life cycle by blocking its interaction with the host [10,11]. An increasing number of recent studies have used computational methods for identifying new drug targets or drug repurposing candidates. Nowadays, various natural compounds have been repurposed/tested using computer-aided drug discovery programs [12,13]. In silico approaches are cost-efective and actionable approaches that can be used to screen several compounds, enabling the discovery of drug combinations and novel drugs.
Flavonoids are natural phytochemicals with antiviral properties which have been discovered to inhibit diferent targets of SARS-CoV [14] such as interfering with spike protein and blocking enzymatic activities of viral proteases, 3-chymotrypsin-like proteases (3CL pro ), and papain-like proteases (PL pro ). Previous research suggests that methylated favonoids, such as retusin, could be used as an antiviral or adjuvant medication in the treatment of COVID-19 [15]. Flavonoids are promising plant-derived chemicals for treating SARS-CoV-2 infection, either through direct antiviral efects or by controlling the host immunological response to viral infection [16]. Previously, several research studies have investigated the use of favonoids against SARS-CoV-2 [16]. An in silico technique was used to examine and compare several favonoids known to have anti-infammatory and antiviral characteristics in an attempt to suppress the spike glycoprotein of SARS-CoV-2, revealing naringin as a potential therapeutic candidate for COVID-19 [17]. Based on molecular docking and ADMET analysis, favonoids, such as cyanidin-3-(p-coumaroyl)-rutinoside-5-glucoside, delphinidin-3-O-beta-D-glucoside 5-O-(6-coumaroyl-beta-D-glucoside), albireo delphine, apigenin 7-(6''-malonylglucoside), and (-) maackiain-3-O-glucosyl-6''-O-malonate were demonstrated as potent inhibitors against the spike protein, 3CLpro, and RdRP of SARS-CoV-2 [18]. Furthermore, molecular docking and simulation studies on favonoid-protein complexes, including luteolin-spike protein and mundulinol-spike protein, have been found to demonstrate strong interactions compared to recently used FDA-approved drugs such as favipiravir, hydroxychloroquine, and lopinavir [19]. Flavonoids, such as dorsilurin E, euchrenone a11, and sanggenol O, were found to inhibit SARS-CoV-2 main protease (M pro ) [20]. Similarly, cyanidin and quercetin were found to inhibit the RNA polymerase function and block interaction sites of the spike protein [21,22]. Te 3CL-protease inhibition assay, cytotoxicity study, total favonoid assay, and molecular docking study were all used to examine a mixture of 11 favonol glycosides prepared from S. persica. It was discovered that these compounds inhibited M pro , spike protein fusing onto its allosteric site [23]. In vitro assay of the extract prepared from muscadine grapes showed inhibitory activity against the M pro of SARS-CoV-2 [24]. Te foundation for future research has been laid by earlier studies of favonoids' efects on SARS-CoV and SARS-CoV-2 inhibition.
Flavonoids, an important class of natural products, can afect CoVs at the initial phase of the entrance, replication, and virion release from the host cells [25]. A stable complex between the spike protein and favonoids is expected to form due to ortho di-OH hydroxyl groups in the B ring of favonols [26]. Many favonoids have been found to block the life cycle of multiple CoV targets (spike protein, proteases, TMPRSS2, etc.) through various mechanisms [27,28]. Flavonoids, therefore, can act as prophylactic, therapeutic, or indirect inhibitors [29]. Flavonoids in this research were chosen based on previously reported in vitro assay data (Table S1) to assemble cyanidin and other favonoids as promising drug candidates against SARS-CoV-2. Furthermore, in vitro assay on the crude extract of red grapes has been carried out for the validation of computational predictions.

Preparation and Molecular Modelling of the Target
Protein. Te crystal structure of the SARS-CoV-2 RBD of spike glycoprotein (PDB ID: 7NX8) with a resolution of 1.95Å was retrieved from the Protein Data Bank (PDB) ( Table S2). Te obtained protein was mutated by replacing threonine (T) with asparagine (N) at position 417 through a sequence editor utilizing the protein builder in Molecular Operating Environment (MOE) software (version 2020.0901). Te structure was optimized using the MOE protein preparation module by removing water molecules, adding hydrogen atoms, and assigning atomic charges to all protein atoms. Moreover, energy minimization was done using diferent force feld parameters of the MOE preparation module.

Preparation of Ligands.
After a comprehensive literature review, the favonoids were selected based on their antiviral properties against SARS-CoV and SARS-CoV-2, and the structures were accessed from the PubChem database [30] and Chem Spider [31]. Finally, the chosen favonoids were processed into mol2 fle format using open babel, and the energy minimization for molecular docking was done using the MOE software [32]. Te structures of the selected 36 natural antiviral favonoids and their IC 50 values against several viruses, including SARS-CoV, SARS-CoV-2, and dengue, are shown in Figure 1.

Active Site Prediction. MOE Site
Finder, which is based on the concept of alpha spheres, was used to determine the amino acid residues involved in the active pocket of the spike protein, where the relative positions and accessibility of the receptor atoms were considered [33]. Furthermore, using the site fnder, the number and name of the residues of the active site were predicted. Table 1 shows the detected cavities in the spike protein region. Figure 2 shows the binding cavity of the spike protein, showing the binding site along with binding residues.

Molecular Docking and Validation. MOE and GOLD
(genetic optimization for ligand docking) (version 4.0.1) software packages [38] were used to predict the proteinligand interactions. MOE docking suite uses the virtual screening protocol of MOE in combination with the triangle matcher docking algorithm and London dG scoring function, and GOLD is based on a genetic algorithm. Te docking analysis was performed on a Microsoft Windows workstation (Intel Core i5-9400 CPU processor and system memory 4 GB RAM). A molecular docking study was conducted on selected favonoids with SARS-CoV-2 spike protein that carries the K417N mutation. Te docking results were validated by extracting the commercial drug and topscored metabolites from their original binding site and redocking them into the same position using the GOLD default docking protocol [39]. MOE provides S-score as its scoring function. Te best pose obtained from MOE is further processed into GOLD software, and fnally, the GOLD score is obtained. To validate the docking results, the lowest energy pose obtained on redocking and the previous docking positions of the metabolites were superimposed, and its root means square deviation (RMSD) was calculated.

Molecular Dynamics Simulation.
Te GROMACS tool was used to produce the protein topology fle and parameterize the natural favonoids' ligand topology. MD simulation for 100 ns was performed in GROMACS 5.1.1 using GROMOS 43a1 force feld for the protein-ligand system and HSA-ligand to clarify the facts behind the efciency of this ligand in protein inhibition [40]. Te PRODRG server was used to obtain a required fle of ligand [41]. During the initial step of the simulation, a cubic box was generated around the protein-ligand complex and solvated with Simple Point Charge (SPC) at the range of 1.0 nm between the wall of components of the protein complex [42]. Sodium ion (Na + ) and chloride ion (Cl − ) with an ionic strength of 0.1 M were added to neutralize the solvated system. Te system's energy was minimized through 50,000 steps of the steepest descent approach, followed by an equilibrium process using Berendsen thermostat (NVT) ensembles. Te integral time phase was two femtoseconds (fs), and the neighbor list was updated for every twentieth step with a cut-of range of 12Å using the grid option. Periodic boundary conditions (PBCs) were used with a constant number of particles in the system, constant pressure, and constant temperature (NPT) simulation criteria. Equilibration of the system at 1 bar pressure for 1 ns was connected using a Parrinello−Rahman barostat in this simulation [42,43]. Using trajectories obtained from MD simulations, root mean square Deviation (RMSD), root mean square fuctuation (RMSF), and the radius of gyration (Rg) were analyzed. Te average structure of the complex was estimated within the last 10 ns trajectory of MD simulations before the RMSF computation, and then each residue around the ligand was aligned to the average structure. During the simulation, the stability of the C-alpha atoms in amino acids was considered. Te stability of the MD trajectories was investigated using the backbone RMSD values of atoms about the spike protein complex, and the time evolution plot of Rg was computed to estimate the conformational stability of the protein-ligand. Te RMSF of all the amino acids around the ligand at 1 nm was determined using the VMD software to investigate the conformational fexibility of the leading active site during the simulation procedure. Te MD runs were executed on an AMD processor (32 cores/64 threads) with 128 GB RAM. Te visual analysis of RMSD, RMSF, and Rg was done using VMD.

Calculation of Binding Free Energy (BFE).
Te prime MM-GBSA [44] of the Schrodinger suite with the OPLS force feld is a popular method to calculate the BFE of binding ligands to proteins. Tis approach is based on docking a ligand and a protein and calculating binding energy using the following equation [45]: where ∆G bind is the BFE of the protein-ligand system and ∆E is the diference in the minimized energies between the protein-ligand complex and the sum of the energies of the free protein and ligand. ∆G Solv is the diference in the GBSA solvation energy of the protein-ligand complex and the sum of the solvation energies for the free protein and free ligand. ∆G SA is the diference in surface area energies for the complex and the sum of the surface area energies for the free protein and free ligand. To prioritize the lead inhibitors, the MM-GBSA technique was applied as a rescoring function.
To optimize the molecules and choose the best compounds, BFE obtained from prime MM/GBSA calculations was considered, along with the docking scores.
2.8. In Vitro Spike S1-RBD and ACE2 Inhibitory Activity of SARS-CoV-2 by Enzyme-Linked Immune Sorbent Assay (ELISA). Te crude methanolic extract of red grape (Vitis vinifera) was prepared with 32.5% yield by the maceration method described elsewhere [46]. A 96-well plate coated with recombinant 2019-nCoV S1-RBD [47] (catalog no. PRncov-2-PL, NOVATEINBIO) was added with increasing concentrations of crude extract to determine whether it inhibits the interaction between hACE2 and SARS-CoV-2 S1-RBD. For each concentration, the assay was performed in triplicate, and % inhibition was expressed as a mean-± standard error of the mean of triplicate. After incubation for 2 hours at 37°C and then washing for 3×PBS (pH 7.2), plates were blocked with 1% BSA and 0.05% Tween-20 in PBS. After washing 3 times, 100 μL hACE2 receptor protein (catalog no.: PR-nCoV-4) (0.1-0.2 μg/mL) was added and incubated for 1 hour in the binding bufer (0.1% BSA in PBS, pH 7.2). Ten, 100 μL goat anti-human IgG-Fc, HRP conjugated 1 : 500 in binding bufer was added after washing plates three times. Again, after three washes by washing bufer, 3,3′,5,5′-tetramethylbenzidine (TMB) was added for the signal and the reaction was stopped by adding an acidic solution. Te plates were read at 450 nm for absorbance, and % inhibition was calculated using GraphPad.

Pharmacokinetic Profles of Flavonoids.
All 36 favonoids were chosen to have a high absorption rate. In addition, their skin permeability, the volume of distribution at steady-state (VDss), CNS permeability, and blood-brain barrier (BBB) permeability were investigated because they play a crucial role in determining drug distribution. Among diferent cytochromes P450 (CYPs) enzymes, the main focus of our study was human cytochrome P450 3A4 (CYP3A4), which was found to be inhibited by favonoids 1, 2, 4, 6, 10, 11, 16, 19, 23, 26, 28, 32, 33, 34, and 35 indicating that they may be metabolized in the liver shown in Table S3. Te toxicity of the selected favonoids was also predicted using ProTox-II, which computes median lethal dose (LD 50 ) values and toxicity classes. . In this investigation, most of the favonoids passed Lipinski's criteria, indicating that they are safe to use as drugs. Table S5 shows the ADME molecular descriptors of the selected favonoids designed to inhibit SARS-CoV-2 by the SwissADME server. Table 2 shows the fve key physiochemical parameters of the selected 6 favonoids which are used to test Lipinski's rule of 5 to evaluate drug-likeness, and Table 3 shows the ADMETprofles of the selected 6 favonoids.

Molecular Docking Analysis.
As indicated earlier, we used the K417N mutant of SARS-CoV-2 spike protein in molecular docking analysis as the mutation at the 417 position by asparagine (N) reduces the therapeutic efcacy of a combination of monoclonal antibodies [6]. Te IC 50 values, GOLD ftness score, bond length, and spike protein interacting residues for the top-scored favonoids along with FDA-approved drugs, lopinavir, and remdesivir, are presented in Table 4. A prior in vitro investigation on SARS-CoV-2 identifed cyanidin as a promising treatment candidate, and the current study shows that it has an appropriate GOLD ftness score of 51.91, indicating its potency as the spike protein inhibitor. Furthermore, 4′-O-methyldiplacol, mimulone, neobavaisofavone, malvidin, and tomentin E also interact appropriately with the binding site of the spike protein with GOLD ftness scores of 63.83, 61.60, 53.90, 52.01, and 50.91, respectively. As shown in Figure 3, cyanidin interacts with the spike protein through Asn 343, Ser 371, Asn 437, Asn 440, and Ser 372 via hydrogen bonds and hydrophobic interactions. Similarly, 4′-O-methyldiplacol, mimulone, malvidin, tomentin E, and neobavaisofavone interact with the protein through diferent residues shown in Table 4. Te GOLD ftness score, interacting residues of the protein, and interaction distances of all remaining favonoids studied are displayed in Table S6. Similarly, Figure S1 depicts the 2D and 3D structures of the remaining 4 potential docked protein-ligand complexes.   Figure 4(a)). Te average RMSD fuctuation for the protein and ligand in the spike-cyanidin complex is 0.39 nm, with equilibrium after 80 ns. Te RMSD of cyanidin was found to be stable frst from 36 to 60 ns and later from 80 ns to 100 ns. Similarly, large deviations in RMSD were observed for the HSA-cyanidin complex, indicating an unstable nature of the complex thus formed. Tese observations support that cyanidin binds with the spike protein rather than HSA. Te RMSF was measured to further compute the residual fexibility over 100 ns shown in Figure 4(b). Te RMSF is less than 0.36 nm for each residue surrounding the ligand in the spike protein complex, while for HSA, it is ∼0.5 nm. Te Rg trajectory of spike-cyanidin reaches equilibrium at ̴ 35 ns and is steady during 35-100 ns, indicating that the ligand is efectively ftted at the active site of spike protein [ Figure 4(c)]. In contrast, the Rg trajectory of HSA-cyanidin optimizes equilibrium during ∼60-80 ns. According to Rg plots, the structural compactness of spike protein-ligand remains stable with an average Rg value of 1.6 nm; in contrast, the HSA-cyanidin complex indicates instability with an average Rg value of 3.5 nm. Conclusively, the lower the Rg value, the higher the compactness of protein-ligand, resulting in a stronger interaction between them [53,54].

In Vitro Assays.
Te red grapes crude extract containing cyanidin [55,56] was in vitro tested using enzyme-linked immunosorbent assay in response to promising in silico results that revealed that cyanidin could inhibit the interaction between hACE2 and SARS-CoV-2 S1-RBD. A series of concentrations of extract was incubated with hACE2 on 96well plates coated with S1-RBD, and the signal was generated using an HRP-linked secondary antibody. Notably, we determined that red grapes extract had the potential to inhibit the interaction of hACE2 with S1-RBD protein. Te detail of the measurement of absorbance and percent inhibition for binding of human ACE2 to S1-RBD protein is shown in Table S7. About 50% reduction in the binding of hACE2 with S1-RBD in presence of red grapes crude extract was achieved at approximately 1.25 mg/mL ( Figure S2). Altogether, this approach provides signifcant mediating evidence that cyanidin, a major favonoid in red grapes (Vitis vinifera), may inhibit the binding of hACE2 and SARS-CoV-2 S1-RBD.

Discussion
Despite continued eforts from various research groups, no efective drug against COVID-19 has been developed yet. In the search for a safe and efective drug, the computational approach provides a good platform for the virtual screening of metabolites to select potential drug candidates [57]. Terefore, repurposing natural secondary metabolites could be the most efective way to discover a remedy for COVID-19 infection [58]. To predict the pharmacokinetic properties of potential drugs, ADMET is an important parameter to be considered. For efective metabolism and activity, a good drug candidate should be absorbed in a specifc time frame   [59]. ADMET properties of metabolites are analyzed and short-listing of metabolites is done by using the Lipinski rule of fve, which determines their drug-likeliness property. According to the Lipinski rule, metabolites with a drug-like property should have a molecular weight (mol·wt.) ≤ 500, hydrogen bond acceptor (HBA) ≤ 10, and the hydrogen bond donor (HBD) should be ≤ 5 [60]. Due to the structural diversity of favonoids, they exhibit versatile biological benefts such as anti-infammatory, neuroprotective, and antioxidative, as well as antiviral properties [17]. Diferent specifc favonoids have been reported effective as antiviral agents for some viruses such as dengue, hepatitis C virus, and infuenza by the modifcation of viral proteins, inhibition of viral entry, and inhibition of viral neuraminidase [61]. Te most investigated targets for favonoid inhibition are proteases (3CL pro and PL pro ) of SARS-CoV and MERS-CoV [62]. Flavonoids such as kaempferol are reported to inhibit cation-selective channels formed by the ORF 3a channel of SARS-CoV. As a result, it ultimately inhibits the virus release and stands out as the root of the evolution of new medicinal antiviral drugs [63]. SARS-CoV-2 S-RBD interacts with host receptor ACE-2 via RBD, laying the groundwork for viral entrance in receptor cells [64]. Te emerging variants of SARS-CoV-2 are connected to the mutation of its spike protein, which is important for the entrance of the virus into the host cells [65]. So, our study considered the spike protein as the target protein for molecular docking of favonoids having an Advances in Pharmacological and Pharmaceutical Sciences antiviral activity to analyze the GOLD score and binding interactions of metabolites. Te GOLD score of 4′-Omethyldiplacol, mimulone, and cyanidin was found to be 63.83, 61.60, and 51.91, respectively. Te greater the GOLD score, the greater its protein-ligand interactions [66]. Te strong interaction of selected ligands to the protein can result in the conformational change of the viral protein and hence, inhibit its replication or growth. However, we considered cyanidin as a promising candidate among the previous three, as it follows Lipinski's rule of 5 and has been revealed to inhibit the SARS-CoV-2 protein in its derivative form with an IC 50 value of 65.1 ± 14.6 μM [48]. In addition, cyanidin shows minimum binding energy with the target protein. It is also reported that cyanidin has around 92% inhibition at 150 μM for M pro enzymatic assay [67]. To assess MD trajectories, the RMSD is a crucial calculation. Using the RMSD of a protein-ligand complex system, the average distance generated by the dislocation of a given atom over time is calculated [70]. Te RMSD of cyanidin was determined to be stable frst from 36 to 60 ns and later from 80 ns to 100 ns. Te RMSF is used to calculate the fexibility among amino acid residues. It is important for monitoring local protein alteration since it allows the calculation of the average change detected over a large number of atoms to determine the displacement of a given atom relative to the reference structure [71]. Te binding pocket was found relatively stable during the MD simulations revealed through obtained RMSF for each residue surrounding the ligand in the protein complex. Rg refects the ligand-protein complex compactness, with a smaller radius of gyration indicating a more compact structure [72]. From the calculation of the Rg, cyanidin efectively binds with the spike protein and found that the complex system was suitably stable ranging from 1.675 nm to 1.70 nm. In addition, the MD simulation analysis of the spike proteincyanidin complex was compared with the HSA-cyanidin complex which reveals the stability of the cyanidin complexed with the target protein. It concludes that the cyanidin comes to bind with the spike protein rather than HSA.

Conclusion
We have used a computational chemistry protocol to identify the most promising favonoids with in vitro IC 50 values that inhibit SARS-CoV and SARS-CoV-2 activity. Tis protocol includes ADMET analysis, molecular docking, molecular dynamics simulations, and BFE calculations to predict whether these favonoids are suitable for anti-COVID-19 therapy. Tis research further enhances that cyanidin interacted with SARS-CoV-2 spike protein, exhibiting the best binding poses and forming stable protein-ligand complexes and lower binding-free energy. Moreover, 4′-O-methyldiplacol and mimulone also showed promising data, based on molecular docking analysis, but as it falls under toxicity class IV, this study repelled to take them as potent agents. Additionally, in vitro study on the crude extract of red grapes ensured the candidacy of   Table S1: plant secondary metabolites (favonoids) with their antiviral activity. Table S2: target protein with PDB ID, resolution, and description of the protein selected for docking with complexed inhibitor. Table S3: ADMET properties of selected favonoids using pkCSM webserver. Table S4: prediction of toxicity of favonoids inhibiting metabolic enzymes using ProTox-II. Table S5: ADMET molecular descriptors of selected favonoids designed to inhibit SARS-CoV-2 by swissADME webserver. Table S6: GOLD ftness score, binding energy, and protein-ligand interaction of natural metabolites with spike protein RBD region. Table S7: measurements of absorbance at 450 nm and calculation of % of hACE2 bound to the S1-RBD, detected by an anti-Human HRP antibody and TMB (n = 3). Figure S1: 2D and 3D structures of malvidin, tomentin E, and neobavaisofavone complexed with SARS-CoV-2 spike protein. Figure S2: binding curve of hACE2 receptor to the S1-RBD protein of SARS-CoV-2 in the presence of a wide concentration range of crude extracts from red grape as determined by ELISA. (Supplementary Materials)