Phytochemicals as potential inhibitors for COVID-19 revealed by molecular docking, molecular dynamic simulation and DFT studies

The COVID-19 pandemic outbreak demands the designing of potential drugs as there is no specific treatment available. Thanks to their safety and effectiveness, phytochemicals have been used to treat various diseases, including antiviral therapeutics. Molecular docking is a simple, quick, and effective way to screen a variety of molecules for structure-based drug design. Here, we investigate molecular docking experiments on compounds present in plant species, Cocculus hirsutus and Rhodiola rosea and show their potential for the treatment of COVID-19. Almost all the components showed higher binding affinity than the built-in ligand, and those with significantly higher binding affinity were explored further. Molecular mechanics-based generalized born surface area calculations were used to re-rank the top candidates, rhodionidin and cocsoline, and their stability in association with viral protease was confirmed. Density functional theory was used for detailed investigations of the geometries, and electrical properties of rhodionidin and cocsoline. Using the frontier molecular orbitals method, the charge transfer within the molecule was calculated. Chemical reactivity and intermolecular interactions were studied using molecular electrostatic potential maps. These in silico discoveries will simulate the identification of powerful COVID-19 inhibitors, and similar research is likely to make a significant contribution to antiviral drug discovery. Supplementary information The online version contains supplementary material available at 10.1007/s11224-022-01982-4.


Introduction
Infectious diseases pose major threat to public health and social well-being. The coronavirus family gained immense attention, when a new human coronavirus called COVID-19 first appeared in the city of Wuhan, China in December 2019 [1]. The global epidemic of severe acute respiratory syndrome (SARS-CoV-2) has impacted the whole world economy as well [2]. Coronavirus comes under the Coronaviridae [3] family in the order nidovirale [4]. Coronavirus contains crown-like spikes on the virus' outer surface [5]. Alpha, Beta, Gamma, and Delta coronaviruses are genetically evolved strains of the reference 2019 SARS-CoV-2 genome and are classified as Variants by WHO [6]. The nucleotide sequence present in SARS-CoV-2 indicates that it belongs to the betacoronavirus category [7,8]. Symptoms of COVID-19 include dry cough, fever, malaise, shortness of breath, and breathing difficullty [9,10]. Drugs such as hydroxychloroquine [11], choloroquine phosphate [12], lopinavir [13], and remdesivir [14] are the antiviral therapies available for the coronavirus as of now. However, their effectiveness remains uncertain and needs further evaluation. Several therapeutic advancements are underway in order to find an effective medication that can directly target viral proteins.
Plant species are immense source of biologically important secondary metabolites such as flavonoids, alkaloids, glycosides, tannins, terpenoids, polyphenols, and saponins [15,16]. Considering the multiple biological activity, researchers are focusing on the several plant-based compounds that can act as potent drug components for various diseases including SARS-CoV-2. However, detection of potential inhibitors using traditional methods is expensive as well as time-consuming. Therefore, in recent years, the use of in silico identification of druggable molecules is gaining tremendous attention as they offer fast and easy detection of potent drugs against various diseases [17,18]. The discovery and isolation of medicinally important therapeutic compounds from higher plants for various diseases have been well established [19][20][21][22][23].
The screening of herbs that may contain anti-coronavirus compounds from a huge library of species is a big challenge. Cocculus hirsutus L.(C.hirsutus), a climbing shrub abundant with alkaloids, flavanoids, and phenolic compounds, comes under Menispermaceae family [24]. The plant has historically been disparaged for its peculiar property of healing many kinds of cuts and wounds and is common in tropical and subtropical environment [25]. It has been shown that this plant is one of the best candidates for treating various diseases such as gonorrhoea, spermatorrhoea, urinary disorders, diarrhoea, and hyperglycemia [26]. The phytochemicals present in this plant are involved in antimalarial [27], antibacterial [26], hyperglycaemic [28], anti-inflammatory, analgesic [28], diuretic, and laxative [29] activities. The genus Rhodiola, belonging to the family Crassulaceae, is composed of almost 200 species. The well-known species among them is Rhodiola rosea L(R. rosea). This plant is also called golden root or arctic root [30], is mainstream among herbal medicines, and is commonly seen at higher altitudes in the arctic and in mountain ranges throughout Europe and Asia [31]. The medicinal application of R.rosea includes its ability to control psychological stress and mental strength, change the neurotransmitter levels and Central Nervous System (CNS) activity, resistance to high altitude sickness, treat fatigue [32][33][34][35][36], and also act as an anti-depressant and anti-inflammatory drug [37]. Phytochemical analysis on R.rosea showed that the plant is well-equipped with bioactive components like organic acids, flavonoids, monoterpenes, triterpene, and tannins, large amounts of phenolic compounds, and some specific phenylpropane derivatives like rosavins [38].
The key druggable target of SARS-CoV-2 involves 3-chymotrypsin-like protease (3CL pro ) [39,40], papainlike protease (PL pro ) [41][42][43][44], RNA-dependent RNA polymerase [45], and spike (S) proteins [46,47]. In this work, 6LU7, main protease (M pro ) of COVID-19, crystallized by Liu et al. (2020), is used as a potential protein target for the phytochemicals selected. The in silico docking model is the main protease, which regulates the major functions of the virus and has a strongly preserved catalytic domain from the SARS virus [48]. Some of its roles involve the virus' replication, making it an attractive target for drug development [49]. Usable bioinformatics techniques such as molecular docking can be implemented using comprehensive 3D structures of related biomolecule and phytochemicals to rapidly classify promising and potential drug candidates. In order to completely unravel the fundamental interactions responsible for drug-receptor binding and subsequent structural modifications, this can be accompanied with atomistic molecular dynamics, leading to changes in biological system function. The current work focuses on the identification of bioactive components present in the herbs C.hirsutus and R.rosea and their antiviral activity towards COVID 19 main protease 6LU7 with the aid of molecular docking and molecular dynamic simulation studies. Towards the end we use density functional theory (DFT) at the B3LYP/6-31 + + G(d,p) level of theory to find the optimal structures of the molecules from each plant that have the maximum binding affinity with the specified target. Following that, their reactivities were predicted at the same level of theory using frontier orbital investigations. Furthermore, molecular electrostatic potential surfaces were investigated to determine which parts of a molecule are the most reactive nucleophilic and electrophilic against reactive biological potentials.

Preparation of protein/receptor
The crystal structure of COVID-19 main protease in complex with an inhibitor N3 (PDB ID-6LU7) with resolution 2.16 Å was downloaded from RCSB protein data bank database (www. rcsb. org) as.pdb format. The 6LU7 protein contains two chains, A and C, in which A represents the protein and C corresponds to the inbuilt ligand. Prior to the detailed docking analysis, chain A was chosen as the target for this study and then the preparation of the protein was performed using AutoDock tool. The selected protein complex contains an inbuilt ligand N-[ Both the water molecules and non-interacting ions, along with the native ligand, were eliminated from the crystal structure. The missing hydrogens were added in order to alleviate the stress of the crystal structure and make the protein accessible for use in the AutoDock docking simulation program. After the structural minimization, the protein was prepared using AutoDock Tools (version1.5.6) (ADT), graphical user interface, which involves the addition of hydrogen atoms, Gasteiger charge calculation, and merging of non-polar hydrogens to carbon atoms. The generated macromolecular structure was then saved as pdbqt file.

Ligand preparation
Most of the ligands were downloaded from the PubChem database (https:// pubch em. ncbi. nlm. nih. gov/). The downloaded ligands in.sdf format were converted into.pdb file using the Open Babel software (version 2.4.1). The ligand structures that were unavailable in the PubChem database were created using Maestro and the structures were downloaded in pdb format. Ligand input files for docking was prepared using AutoDock Tools and saved as pdbqt files.

Docking
For the detailed analysis of intermolecular interactions and binding modes between the protein and specific phytochemicals present in C. hirsutus and R. rosea, the molecular docking studies were executed for selected ligands using the AutoDock Vina software. A grid box with dimension X:15, Y:22, Z:20 Å with a grid spacing of 1.0 Å centered at X: −9.491, Y:11.784, Z: 65.363A 0 was identified as protein target docking site. The interactions between active sites in the target and ligand conformation, including the type of interaction and bond distances, were identified using Discovery Studio Visualizer.

Binding energy calculations using prime/MMGBSA studies
The calculation of free energy provides a quantitative estimate of protein-ligand interactions that can predict the stability of protein-ligand complexes. Molecular mechanicsgeneralized born surface area (MM-GB/SA) was conducted with Prime module of Schrodinger 2018 to calculate binding energy of protein-ligand complexes. MM-GBSA method in Prime was used for rescoring the docked pose of ligands together with calculation of coulomb energy, covalent binding energy, hydrogen-bonding correction, pi-pi packing correction, lipophilic energy, generalized born electrostatic solvation energy, and Van der Waals energy.

Molecular dynamics simulation
To forecast the stability of protein-ligand complexes, molecular dynamic simulations were performed in Desmond package of the Schrodinger suit. The chosen protein-ligand complexes were first immersed in the TIP4P water box, extending 10 Å beyond any atoms of the complex. The orthorhombic boundary condition box is selected, and buffer method is used for box size calculation. To neutralize charges, counter ions (3 Na + ) were added. The MD was conducted at a temperature of 310 K and 1.63 bar pressure over 100 ns in the NPT ensemble with recording intervals of 12 ps for trajectory. Simulations were performed with the force field of OPLS-3e. The Desmond simulation interaction diagram method from Maestro was used to sketch plots and figures.

Quantum chemical analysis
Density functional theory (DFT) is an intriguing method for probing many biological features of interest. The current study employs DFT analysis to highlight the electronic features of the phytochemicals. The geometry optimizations of the selected compounds cocsoline (from C.hirsutus) and rhodionidin (from R.rosea) were carried out using Gaussian16 at the B3LYP/6-31 + + G(d,p) level. Based on the optimized structures in the gas phase, electronic characteristics such as HOMO-LUMO energies and molecular electrostatic potential (MESP) were computed. The natural population analysis, as well as global reactivity descriptors, was further investigated for the selected compounds.

Results and discussion
This section is organized as follows. First, we identified the chemical components of C.hirsutus and R.rosea and the identified components were docked with 6LU7, the SARS-CoV-2 target protein. The components with highest binding affinity obtained using docking were re-ranked using MMG-BSA calculation. Then, we evaluated the stability of the complex formed between the components having the highest binding affinity with 6LU7, using RMSD values obtained from the MD simulation analysis. The position of binding modes of components, rhodionidin and cocsoline (with the highest binding affinity with 6LU7) was verified using the MD modelling counterpart. Finally, the density functional theory calculations were performed to investigate the structure and reactivity of top ranked candidates.

Molecular docking analysis of phytochemicals present in C.hirsutus and R.rosea
The essential parts of the structure-based drug designing phase involve the docking of small molecules into the receptor's active site and evaluating the binding affinity of ligand towards target protein. In the current study, the selected protein, 6LU7 contains a peptide like inbuilt ligand N3 ( Fig. 1) that exhibits very potent inhibition against COVID-19 virus M pro [50]. In order to compare the binding affinity as well as the inhibition potency of selected phytochemicals, the docking process was also performed on the native ligand.  Table S1 and the binding affinity values of the components with high inhibition efficiency than that of in-built ligand N3 are shown in Fig. 3a. The binding affinity of N3 with the target protein obtained from the docking study was −7.0 kcal/mol. The lowest binding energy value of −9.1 kcal/mol was observed for the phytoconstituent, cocsoline. The flavonoids rutin, liquiritin, and quercetin with binding energy values −8.8, −7.5, and −7.5 kcal/mol respectively prove that they could also act as possible inhibitors against 6LU7 protein target. Isotrilobine, coclaurine, cohirsine, and lupeol, exhibiting binding affinity values like −8.6, −7.6, −7.4, −7.3, and −7.1 kcal/mol respectively, are still higher than that of binding affinity of native ligand N3.
Docking analysis of phytochemicals of the plant C.hirsutus shows that all of them show strong interactions with the surrounding amino acid residues. The two-and three-dimensional interactions of components of C.hirsutus having low binding energy values < −7.4 kcal/mol as well as the inbuilt ligand N3 with 6LU7 are represented in Table 1. Cocsoline, having the highest binding affinity, possess hydrogen bond interaction between its hydroxyl group and amino acid residues     Table S3. The majority of the chemicals found in R.rosea showed good binding affinity values with the specified protein target, 6LU7, according to data from molecular docking. The most of the binding affinity values were found to be greater, with some even being comparable to the native ligand N3. The molecular docking results of bioactive components (along   belonging to the phenyl proponoid group, also showed a high affinity towards the target protein. Docking studies revealed that rosin and vimalin (belonging to phenylproponoids) and rhodioloside (belonging to phenylethanoids) show a binding energy value with the SARS-CoV-2 which is equal to the binding energy of inbuilt ligand with the same target. Other components like picein, rosiridin, gallic acid, rosiridol, geraniol, cinnamyl alcohol, limonene, alphapinene, and decanol exhibit low inhibition for the target. More interestingly, these values are still comparable with the native ligand.
The three-and two-dimensional binding interactions of best scoring components present in R.rosea having low binding energy values < −7.9 kcal/mol with 6LU7 are provided in Table 2. The major interactions involved are conventional hydrogen bonding, alky, pi-alkyl, pi-pi, pi-S, pi-cation, carbon-hydrogen bond, and pi donor hydrogen bond and all of the chemical components shows very strong interaction with the protein studied.
The detailed visualization of stable conformer of rhodionidin, the top-ranked (binding energy value of −9.6 kcal/ mol) SARS-CoV-2 inhibitor at the binding site of 6LU7, reveals the existence of a strong hydrogen bond interaction between hydroxyl groups present in the ligand and GLY143, SER144, LEU141, and HIS163 amino acid residues of the target protein. A π-π interaction exists between the inhibitor and HIS41 residue present in the active site of 6LU7. Alkyl and aromatic fragments present in the inhibitor form hydrophobic interactions such as alkyl and pi alkyl interaction with LEU27, MET49, PRO168, and MET165 residues. Furthermore, the molecule is well fixed within the cavity of the active site.
From the detailed interaction studies of phytochemicals present in the herb R. rosea showing high binding energy with selected protein proved that all components are potential inhibitors and possess strong interactions with the target. The structure, 3D images of the phytochemicals within the active site of 6LU7, and 2D interaction of the phytochemicals present in R. rosea, with binding energy values less than -8.0kcal/mol with the target, are given in Table S4.

Molecular dynamic simulation using Desmond package
Since docking is a static view of the molecule's binding position in the protein's active site, molecular dynamics (MD) simulation attempts to quantify atomic motions over time by integrating the classical motion equation of Newton. The dynamic behavior of the molecular system is simulated in MD, in order to assess the stability of the protein-ligand complex. For the MD analysis with the OPLS3e force field, the docked conformers of cocsoline and rhodionidin with the highest predictive binding energy of −9.1 and −9.6 kcal/ mol were therefore used. The dynamic features of these two selected systems were calculated during the period of 100 ns simulation. We evaluated the stability of the docked compounds in the active pocket of 6LU7 by using root mean square deviation (RMSD) and its effects on the stability of the entire system. RMSD plot in Fig. 5a indicates the rhodionidin:6lu7 complex MD trajectory for 100 ns. The complex appears to be stabilized within a short time with respect to the reference frame at time 0 ns, during the simulation phase. The RMSD variation of 6lu7 was only around 1.6 Å which can be considered insignificant. This in turn means that the protein RMSD appears to be stable during the simulation, indicating that the protein has not undergone major conformation changes. Fluctuation within 20-40 ns of the trajectory was observed in the case of ligand RMSD which indicate major conformational changes of the ligand. After 40 ns, the ligand RMSD stabilized again and remained constant until the end of the simulation period of 100 ns. The ligand has undergone conformational changes with an overall RMSD of 6 Å, possibly because of a large number of rotatable bonds that contribute to its flexibility. The protein RMSD appears to be constant with negligible variation which means that the protein is undergoing minimum conformational change during the simulation. The rhodionidin:6lu7 complex's overall RMSD plot shows that the ligand is less stably connected to the protease binding site.
The RMSD plot in Fig. 5b shows that the cocsoline:6lu7 complex stabilized immediately after the simulation started, relative to the reference frame at 0 ns. During the simulation, the complex appears to be stabilized with regard to the reference frame at time 0 ns. It is possible to see a small divergence at the end of the simulation, i.e., 80 ns. The fluctuation, however, lies within the acceptable range of 1-3 Å, and can therefore be considered non-significant. In this complex, since the ligand and protein backbone RMSD plots were lying over each other, it is possible to conclude a stable complex formation. RMSD value analysis reveals that the complex cocsoline:6lu7 was more stable compared to the rhodionidin:6lu7 complex.
The position of the binding modes of compounds rhodionidin and cocsoline has also been verified by the MD modelling counterpart. During the simulation, protein interactions with the ligand were regulated. Such interactions can be arranged and analyzed as depicted in Fig. 6. Protein-ligand interactions have been classified as hydrogen bonds, hydrophobic, ionic, and water bridges. In order to keep the protein-ligand complex stable, the form of amino acid residues present in the active site of the target protein plays a crucial role. Figure 6a indicates the type of contacts maintained in the entire trajectory of rhodionidin:6LU7 complex. In the case of rhodionidin, GLN189, GLU166, and GLN192 were able to hold down   where one water molecule is used as a bridge between ligand and protein, was also formed with GLU166, GLN189, THR24, THR25, THR26, SER46, ASN142, ARG188, THR190, and GLN192. The various types of interaction between cocsoline and 6LU7 at the active site during the entire simulation period are represented in Fig. 6b. Cocsoline showed hydrogen bonding interaction with THR26 and HIS41 which were 49% and 30% of the time, respectively. It also exhibits weak hydrogen bonding contacts with GLN189, GLY143, HIS164, and GLU166 less than 10% of the time. Cocsoline possesses multiple water bridge hydrogen bonding with THR190, GLN192, GLN189, GLU166, CYS145, GLY143, SER144, TYR54, ASN119, THR25, HIS41, MET49, ASN142, HIS163, HIS164, ASP189, and ARG188 residues.
Relevant timeline contacts of the compounds rhodionidin and cocsoline with the amino acid residues that were present in the targets are also studied and provided in Fig. 7. The darker color reflects the high number of ties to the amino acid.

MMGBSA analysis
MMGBSA is a popular analysis method to calculate the average binding energy of equilibrated MD trajectory of the systems. MMGBSA calculations were conducted in order to provide a better ranking of the ligands, and for the determination of predictive binding energies. The binding free energies of the two selected systems indicate the binding affinity of ligands with 6LU7. The values obtained after measurement were estimated free binding energies with a more negative value, suggesting a stronger binding energy. Table 3 depicts the average binding energy, coulomb energy, covalent binding energy, hydrogen-bonding correction, pi-pi packing correction, lipophilic energy, generalized born electrostatic solvation energy, Van der Waals energy of rhodionidin, and cocsoline with 6LU7. Rhodionidin exhibited the highest negative binding energy of −80.40 kcal/mol. MMGBSA calculation revealed that cocsoline showed an approximate binding energy value of −75.07 kcal/mol.

Geometry optimization
The optimized structures of cocsoline and rhodionidin obtained at B3LYP/6-31 + + G(d,p) level are given in Fig. 8, and the corresponding geometrical parameters are given in Tables S5-S10. The total energy of cocsoline and rhodionidin computed at the same basis set are −1800.92 and −2251.14 Hartree respectively. A molecule's dipole

Natural population analysis and molecular electrostatic potential calculations
Molecular polarizability, dipole moment, electronic structure, molecular reactivity, and many other aspects of molecular systems are influenced by atomic charges. The creation of donor and acceptor pairs involving the charge transferring molecule is suggested by the charge distributions over the atoms. The distribution of electrons in various sub shells of atomic orbitals is revealed by a natural population analysis (NPA) of an organic molecule [52]. The natural charges on atoms of cocsoline and rhodionidin are given in Tables 4 and 5 respectively. The value of natural charges informs that the atom O10 and O5 were showed more electronegativity in cocsoline and rhodionidin respectively than all other oxygen atoms. Among the two nitrogen atoms present in rhodionidin, N7 was more electronegative than N6 atom. On other hand H67 of cocsoline and H73 of rhodionidin were possessed high electropositive value. The three-dimensional charge distributions of molecules will be depicted by the molecular electrostatic potential surface. The molecular electrostatic potential (MEP) surface can also be viewed as a map of electron excess and deficient areas [53,54]. The MEP must be determined to validate the evidence of the drug's reactivity as inhibitors, despite the fact that the MEP indicates the molecule size and shape of the positive, negative, and neutral electrostatic potentials. These could be used to anticipate physicochemical property relationships based on the molecular structure of medications in development. Furthermore, the molecular electrostatic potential can be used to evaluate a drug's responsiveness to electrophilic and nucleophilic assaults. The MEP diagram of both cocsoline and rhodionidin calculated at DFT/ B3LYP/6-31 + + G(d,p) is given in Fig. 9. In the MEP, the largest negative region, which is shown red, is the favored site for electrophilic attack. As a result, an attacking electrophile will be drawn to the negatively charged sites, whereas the blue regions will attract the reverse. The change in the drug's binding affinity with the active site receptor could be due to differences in the mapping of the electrostatic potential around it.

Frontier molecular orbital analysis
The value of highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) energy can be used to determine a molecule's ability to donate and   Fig. 10, the negative phase is represented by the color green, whereas the positive phase is represented by the color red. Rhodionidin showed the most lying HOMO than cocsoline, and consequently, it could be a better electron donor drug. The energy gap of rhodionidin and cocsoline is 3.6662 and 4.7775 eV respectively. Rhodionidin's small energy gap allows electrons to move more freely, making the molecule soft and reactive and this result is in strong agreement with the high dipole moment of rhodionidin. The HOMO of a certain drug and the LUMO with the adjacent residues could share the orbital interactions during the binding process.

Chemical reactivity descriptors
The E HOMO and E LUMO are markers that can be used to forecast a molecule's ionization potential (I = −E HOMO ) and electron affinity (A = −E LUMO ). Other chemical reactivity descriptors such as hardness (η), softness (S), electronegativity (χ), chemical potential (μ), and electrophilicity index (ω) are estimated using the frontier molecular orbitals and are shown in Fig. 11.
The electronegativity (χ = (I + A)1/2) value predicts the molecule's ability to attract electrons, i.e., Lewis acid, whereas lower values of χ indicate a suitable base. The χ value is found to be 3.15 and 4.09 for cocsoline and rhodionidin, respectively. The global hardness (η = (I −A)1/2) is a measure of their charge transfer prohibition, while the global softness (S = 1/2η) characterizes a molecule's ability to accept electrons. Because they may easily transmit electrons to acceptors, soft molecules have a smaller energy gap between frontier molecular orbitals and are more reactive than harder molecules. The hardness of cocsoline (2.39 eV) was higher than that of rhodionidin (1.83 eV) whereas the softness of cocsoline (0.21 eV) is lower than that of rhodionidin (0.27 eV). The electrophilicity index (ω) = μ 2 /2η determined for cocsoline and rhodionidine is 2.07 and 4.56 eV, respectively. The chemical potential (μ = −χ) is found to be −3.15 eV for cocsoline and −4.09 eV for rhodionidin.

Conclusion
Components present in the two medicinal herbs R.rosea and C.hirsutus were subjected to molecular docking studies in order to identify possible COVID-19 M pro inhibitors. Using the AutoDock Vina software, the binding energy values of all the components present in both plants were determined. The findings of in silico molecular docking studies of R.rosea and C.hirsutus indicate that most of the phytochemicals have moderate to potent antiviral activity against SARS-CoV-2 M pro . The components' high affinity is also found to be dependent on the nature and strength of bonding within the protein's active site. The flavonoid rhodionidin, which is found in the R.rosea plant and has the lowest binding energy of −9.6 kcal/mol, and the alkaloid cocsoline, which is found in the species C.hirsutus and has a binding energy of −9.1 kcal/mol, have been described as the best potential drug against the novel coronaviruses. MMGBSA analysis was used to help rerank these ligands showing the highest binding affinity for 6LU7. MD simulations of rhodionidin and cocsoline with the target indicate that the complex cocsoline:6lu7 was more stable compared to the rhodionidin:6lu7 complex. The DFT/ B3LYP/6-31 + + G(d,p) approach was used to optimize the molecular structures of cocsoline and rhodionidin (having highest binding affinity with the target protein), as well as to identify their geometrical characteristics. Frontier orbitals, gap energies, natural population analysis, and reactivity descriptors are some of the molecular features that have been discussed. According to the findings, the minimal energy gap is associated with a high binding affinity for rhodionidin. The molecule gets softer and more reactive as the gap energy decreases, making electron flow easier. Following that, the estimated MEP maps reveal that positive potential sites are more advantageous for nucleophilic attack, while negative potential sites are more favorable for electrophilic attack. These findings further motivate in vitro and in vivo studies of the identified components for the treatment of COVID-19. These findings are also expected to stimulate further research into the production of safe and successful anticoronavirus or other antiviral drugs derived from naturally occurring compounds.