Conducting the RBD of SARS-CoV-2 Omicron Variant with Phytoconstituents from Euphorbia dendroides to Repudiate the Binding of Spike Glycoprotein Using Computational Molecular Search and Simulation Approach

(1) Background: Natural constituents are still a preferred route for counteracting the outbreak of COVID-19. Essentially, flavonoids have been found to be among the most promising molecules identified as coronavirus inhibitors. Recently, a new SARS-CoV-2 B.1.1.529 variant has spread in many countries, which has raised awareness of the role of natural constituents in attempts to contribute to therapeutic protocols. (2) Methods: Using various chromatographic techniques, triterpenes (1–7), phenolics (8–11), and flavonoids (12–17) were isolated from Euphorbia dendroides and computationally screened against the receptor-binding domain (RBD) of the SARS-CoV-2 Omicron variant. As a first step, molecular docking calculations were performed for all investigated compounds. Promising compounds were subjected to molecular dynamics simulations (MD) for 200 ns, in addition to molecular mechanics Poisson–Boltzmann surface area calculations (MM/PBSA) to determine binding energy. (3) Results: MM/PBSA binding energy calculations showed that compound 14 (quercetin-3-O-β-D-glucuronopyranoside) and compound 15 (quercetin-3-O-glucuronide 6″-O-methyl ester) exhibited strong inhibition of Omicron, with ΔGbinding of −41.0 and −32.4 kcal/mol, respectively. Finally, drug likeness evaluations based on Lipinski’s rule of five also showed that the discovered compounds exhibited good oral bioavailability. (4) Conclusions: It is foreseeable that these results provide a novel intellectual contribution in light of the decreasing prevalence of SARS-CoV-2 B.1.1.529 and could be a good addition to the therapeutic protocol.


Introduction
In October 2021, a new SARS-CoV-2 B.1.1.529 variant emerged in South Africa, designated Omicron by the World Health Organization (WHO). By early 2022, over 100,000 Omicron genomes had evolved as Omicron begun to dominate SARS-CoV-2 infections around the world [1]. Because it contains more mutations than any other strain, it is more transmissible than previous strains. Many of the changes are found in the spike protein, which is involved in the transmission of the virus.
Natural ingredients offer a wide range of chemical properties, including antiviral activity, and thus could be used to treat coronavirus infections. Plants and their secondary metabolites that act against targets associated with SARS-CoV-2 infection could be useful leads for developing drugs against the newly emerged Omicron. The discovery of antiviral drugs and effective therapeutic techniques is lengthy and laborious. Therefore, natural chemicals are often considered as attractive alternative treatment solutions because they are also the main sources of antibacterial and antiviral drugs. Recently, it was suggested that dietary flavonoids could regulate the severity of SARS-CoV-2 disease by affecting ACE2 prevalence and function [2,3]. Several publications have also stated that polyphenols target the renin-angiotensin system by modulating angiotensin II levels in mice [4,5].
In silico studies recently reported that flavonoids such as quercetin and rutin and polyphenols such as epigallocatechin gallate, myricetin, and quercetagetin showed a high rate of inhibition against the RNA-dependent RNA polymerase (RdRp) of SARS-CoV-2 [6].
Euphorbia dendroides L. (family Euphorbiaceae), one of the species with high polyphenol content, is a perennial small tree growing in Sollum and Mersa Matruh in Egypt [7]. The genus Euphorbia is unique in that it includes highly reputed species used in traditional medicine against various human diseases, such as respiratory diseases, inflammation, skin diseases, diarrhea, migraine, gonorrhea, warts, and intestinal parasites and exerts a laxative effect [8][9][10]. In addition to E. dendroides, other plants in the genus Euphorbia and their constituents have also recently gained medicinal importance and are used for various diseases, including as anticancer [11][12][13], antioxidants [14], and antiviral agents, and to target multidrug resistance [15,16] and COVID-19 [17][18][19]. Our previous phytochemical studies on E. dendroides' aerial parts revealed the presence of ten phenolic compounds, including six flavonoids, one phenolaldehyde, and three phenolic acids [20], as well as six cycloartane triterpenes and lupeol triterpenes, in addition to β-sitosterol and three fatty acids [13]. The constituents of medicinal plants may serve as targets for the development of therapeutic candidates against some SARS-CoV-2 proteins [21,22]. Consequently, in this study, seventeen phytoconstituents from the categories of flavonoids, phenols, and triterpenes that we had previously obtained from this plant were evaluated by in silico studies as potential candidates against the Omicron variant of severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2).
In this study, we performed computational molecular screening against the Omicron receptor-binding domain (O-RBD) to explore and develop drugs from natural plant constituents that are effective against the current pandemic virus.

Molecular Docking
An important aspect of the drug discovery strategy is molecular docking analysis, which can be used to specify protein-ligand interactions in the active site of the target protein. In this context, all seventeen compounds were subjected to molecular docking calculations to investigate their potency as anti-Omicron drugs. For each compound, only the docking pose with the highest docking score was selected from the nine docking poses. The estimated docking scores for all seventeen compounds are shown in Table S2. As shown in the data deposited in Table S2, eight compounds had docking scores of less than −7.0 kcal/mol. For the remaining compounds, three compounds were in the range of −7.2 to −7.7 kcal/mol and six compounds were in the range of −7.9 to −8.8 kcal/mol. The average docking score for the seventeen compounds was calculated to be −7.04 kcal/mol.
Remdesivir is a nucleotide analog prodrug and has been recently subjected to in vitro experiments as an anti-Omicron drug [34]. In order to assess the potentiality of the discovered compounds, the binding features and affinities of the top six compounds were compared to those of remdesivir against Omicron. Table 1  To recognize the other types of interactions, two-dimensional representations of the interactions of these six potent compounds with the main active site residues are displayed in Figure 1. It is worth mentioning that compound 15 exhibited a pi-alkyl interaction with LEU455 (4.96 Å) and ARG493 (4.17 Å). Moreover, it also formed a pi-cation interaction with ARG403 (4.92 Å), as well pi-pi stacked and pi-pi T-shaped with TYR501 (4.24 Å) and HIS505 (4.92, 4.99 Å), respectively. In the case of compound 14, it formed two extra types of interactions, pi-pi T-shaped with HIS505 (4.92, 5.05 Å) and a pi-cation interaction with ARG403 (4.90 Å). Similar to compound 15, compound 16 formed pi-alkyl, pi-pi stacked, and pi-cation interactions with the same residues but different bond lengths. Excluding pi-alkyl interaction, compound 17 showed the same types of interactions that existed in compound 16. Notably, compound 13 was the only one that exhibited no pi-cation interaction. Compound 12, the lowest-ranking of the six compounds in terms of docking score, showed only two additional interactions other than hydrogen bonding. It is worth highlighting that residues ARG403, TYR501, and HIS505 favored the formation of pi-cation, pi-pi stacked, and pi-pi T-shaped interactions, respectively, with RBD residues.

Molecular Dynamics (MD) Simulations
The main purpose of applying molecular dynamics (MD) simulations is to investigate the conformational flexibilities and stabilities of studied protein-ligand complexes. Accordingly, MD simulations up to 200 ns were performed for the best-ranking six complexes. Binding energy calculations were also run using the MM/PBSA approach. Data of estimated binding energies are represented in Figure 2.
Notably, compounds 14 and 13 showed similar binding energies over 50 ns with values equal to −44.8 and −46.4 kcal/mol, respectively. Furthermore, compound 15 and compound 16 showed identical values of binding energy. Interestingly, compound 12, the lowest of the six compounds in terms of docking score, also exhibited the lowest binding energy, with a value equal to −18.6 kcal/mol.      pi-pi stacked, and pi-cation interactions with the same residues but different bond lengths. Excluding pi-alkyl interaction, compound 17 showed the same types of interactions that existed in compound 16. Notably, compound 13 was the only one that exhibited no pi-cation interaction. Compound 12, the lowest-ranking of the six compounds in terms of docking score, showed only two additional interactions other than hydrogen bonding. It is worth highlighting that residues ARG403, TYR501, and HIS505 favored the formation of pi-cation, pi-pi stacked, and pi-pi T-shaped interactions, respectively, with RBD residues.

Molecular Dynamics (MD) Simulations
The main purpose of applying molecular dynamics (MD) simulations is to investigate the conformational flexibilities and stabilities of studied protein-ligand complexes. Accordingly, MD simulations up to 200 ns were performed for the best-ranking six complexes. Binding energy calculations were also run using the MM/PBSA approach. Data of estimated binding energies are represented in Figure 2.

Post-MD Analyses
The structural stability and conformational variations of protein-ligand docked complexes over the MD simulation process can be evaluated by root-mean-square deviation (RMSD) measurement. Lower values of RMSD give an indication of tight binding. The following equation describes how RMSD was estimated.

Post-MD Analyses
The structural stability and conformational variations of protein-ligand docked complexes over the MD simulation process can be evaluated by root-mean-square deviation (RMSD) measurement. Lower values of RMSD give an indication of tight binding. The following equation describes how RMSD was estimated.
where R i is the vector connecting the positions of atom i [of N atoms] in the reference snapshot and the current snapshot after optimal superposition. Obtained results of RMSD analysis over 200 ns MD simulation are plotted in Figure 3. It is worth mentioning that compound 14 exhibited overall stability compared to the other five compounds, with an average RMSD value of 2.84 Å. Comparatively, the RMSD values of the other compounds were relatively high. These results showed that compound 14 was tightly bound and had no effect on the overall topology of Omicron.
where Ri is the vector connecting the positions of atom i [of N atoms] in the reference snapshot and the current snapshot after optimal superposition. Obtained results of RMSD analysis over 200 ns MD simulation are plotted in Figure  3. It is worth mentioning that compound 14 exhibited overall stability compared to the other five compounds, with an average RMSD value of 2.84 Å. Comparatively, the RMSD values of the other compounds were relatively high. These results showed that compound 14 was tightly bound and had no effect on the overall topology of Omicron. Root-mean-square fluctuation (RMSF) analysis was conducted to investigate the flexibility of the Omicron residues over the 200 ns MD simulations. Briefly, higher RMSF values indicate greater flexibility of protein residues, whereas low RMSF values imply limitations of residues' movement, and accordingly less flexibility. RMSF can be estimated using the following equation: Root-mean-square fluctuation (RMSF) analysis was conducted to investigate the flexibility of the Omicron residues over the 200 ns MD simulations. Briefly, higher RMSF values indicate greater flexibility of protein residues, whereas low RMSF values imply limitations of residues' movement, and accordingly less flexibility. RMSF can be estimated using the following equation: The RMSF of the atom i with j from 1 to 3 for the x, y, and z coordinate of the position vector P of the atom and k over the set of N evaluated snapshots was calculated.
The RMSF data are shown in Figure 4.  In order to assess whether the protein-ligand complexes were stably folded or not, radius of gyration (Rg) analysis was performed over the 200 ns MD simulations. Values of Rg give an indication of the compactness of the protein structure within the system [35] A more compact protein structure can be observed in the case of lower values of Rg. During the 200 ns simulation course, all compounds exhibited acceptable behavior of Rg, which can be observed in Figure 5. The range of Rg values for the six compounds was from 18.68 to 19.14 Å. What also can be noticed from Figure 5 is that compounds 14 and 15 showed constant Rg behavior during the simulation, which indicates high compactness of the protein structure.  In order to assess whether the protein-ligand complexes were stably folded or not, radius of gyration (Rg) analysis was performed over the 200 ns MD simulations. Values of Rg give an indication of the compactness of the protein structure within the system [35] A more compact protein structure can be observed in the case of lower values of Rg. During the 200 ns simulation course, all compounds exhibited acceptable behavior of Rg, which can be observed in Figure 5. The range of Rg values for the six compounds was from 18.68 to 19.14 Å. What also can be noticed from Figure 5 is that compounds 14 and 15 showed constant Rg behavior during the simulation, which indicates high compactness of the protein structure.
Solvent-accessible surface area analysis (SASA) was performed to represent the area of protein exposed to solvent. The averaged SASA values of the six structures were 10,501, 10,428, 106,30, 11,007, 10,553, and 10,886 for compounds 12 to 17, respectively. Despite the compact folding of 7QNW, the increase in SASA in compound 15 compared to the other compounds indicated obvious conformational changes due to ligand binding. Consequently, SASA may provide details about the protein's ability to interact. As shown in the data in Figure 6, all six compounds exhibited relative stability and compactness over the course of the 200 ns MD simulation.

In Silico Drug Likeness
Based on Lipinski's rules, physicochemical properties were investigated to understand the studied compounds' molecular features better. The Molinspiration tool https://www.molinspiration.com (30 April 2022) was utilized to compute the in silico molecular features of compounds. The investigated Lipinski's parameters, topological polar surface area (TPSA), as well percentage of absorption (% ABS) of Lipinski's parameters were anticipated and are presented in Table 2.
In order to assess whether the protein-ligand complexes were stably folded or not, radius of gyration (Rg) analysis was performed over the 200 ns MD simulations. Values of Rg give an indication of the compactness of the protein structure within the system [35] A more compact protein structure can be observed in the case of lower values of Rg. During the 200 ns simulation course, all compounds exhibited acceptable behavior of Rg, which can be observed in Figure 5. The range of Rg values for the six compounds was from 18.68 to 19.14 Å. What also can be noticed from Figure 5 is that compounds 14 and 15 showed constant Rg behavior during the simulation, which indicates high compactness of the protein structure. Solvent-accessible surface area analysis (SASA) was performed to represent the area of protein exposed to solvent. The averaged SASA values of the six structures were 10,501, 10,428, 106,30, 11,007, 10,553, and 10,886 for compounds 12 to 17, respectively. Despite the compact folding of 7QNW, the increase in SASA in compound 15 compared to the other compounds indicated obvious conformational changes due to ligand binding. Consequently, SASA may provide details about the protein's ability to interact. As shown in the data in Figure 6, all six compounds exhibited relative stability and compactness over the course of the 200 ns MD simulation.

In Silico Drug Likeness
Based on Lipinski's rules, physicochemical properties were investigated to understand the studied compounds' molecular features better. The Molinspiration tool https://www.molinspiration.com (30 April 2022) was utilized to compute the in silico molecular features of compounds. The investigated Lipinski's parameters, topological polar surface area (TPSA), as well percentage of absorption (% ABS) of Lipinski's parameters were anticipated and are presented in Table 2. As revealed from the data in Table 2, the investigated two compounds showed lower  As revealed from the data in Table 2, the investigated two compounds showed lower values of miLogP, indicating that these compounds possess adequate permeability via the cell membrane. The molecular weights of the inspected compounds did not exceed 500 (calc. 478.4 and 492.4). Hydrogen bond donors (nON) were found to total 13. Furthermore, the number of hydrogen bond acceptors (nOHNH) was found to range between 7 and 8. Although the values of nON and nOHNH were higher than the ideal values, it was reported that this defect did not exert a remarkable effect on the compound's diffusion and transportation, as many FDA-approved drugs transcend the optimum Lipinski values of nON and nOHNH [36]. The estimated %ABS values were in the range of 30.0-35.0%. The TPSA values were also detected in the range of 215.0-230.0, which indicates the high bioavailability of these discovered compounds.

Plant Material
The aerial parts of E. dendroides were collected in April 2017 in Mersa Matruh on the northwestern coast of Egypt. The plant material was authenticated by Dr. Omran Ghaly, with a PhD in Plant Taxonomy, at the Desert Research Center. A voucher sample (CAIH-30-12-2017-R) was deposited in the herbarium of the Desert Research Center, Cairo, Egypt.

Protein Preparation
The three-dimensional crystal structure of the receptor-binding domain (RBD) of the SARS-CoV-2 Omicron variant (PDB ID: 7QNW, resolution: 2.40 Å) was retrieved and used for all in silico analyses. The downloaded viral target was prepared by removing ions, water molecules, and hetero-atoms. To identify the protonation states of the protein residues, the H++ web server [37,38] was used. Accordingly, all missing hydrogen atoms were successfully added. For the H++ calculations, physiological parameters such as pH, salinity, internal dielectricity, and external dielectricity were set to 6.5, 0.15, 10, and 80, respectively.

Inhibitor Preparation
Chem3D Pro 12.0 software (version 12.0.2) was utilized to sketch and analyze the seventeen extracted compounds' chemical structures. All studied compounds were subjected to energy minimization using the MM2 force field. Before a molecular docking study, such an energy minimization step is required to reduce the influence of any potential unfavorable torsion angles, bond angles, bond lengths, or undesirable non-bonded interactions [39]. The names of the investigated compounds and their 2D chemical structures are illustrated in Table S1.

Molecular Docking
Molecular docking is considered the best tool in computational drug discovery to determine the efficacy of the compounds under study [40]. AutoDock Vina was used to study the binding affinities for these compounds [41]. All parameters were left in their default modes in this study, except for the exhaustiveness parameter, which was set to 200. Residues of the O-RBD were enclosed by a docking grid box with XYZ dimensions of 25 × 25 × 25 (Å). In addition, the grid spacing was set to 1.0 Å. The generated nine poses of the docked inhibitors were evaluated and the best one was selected. BIOVIA Discovery Studio was used to visualize the protein-ligand interactions [42].

Molecular Dynamics Simulations
Molecular dynamics (MD) simulations were conducted using the YASARA Structure (version 21.12.19) protocol [43] to obtain a better understanding of the stability of the protein-ligand complexes. Within the MD simulations, the AMBER14 force field was the utilized force field. Execution of the initial energy minimization was performed using the steepest descent algorithm. The MD simulations were conducted for amino acid residues at the default physiological value of pH (7.4). Water molecules were successfully introduced into the system at constant temperature and pressure conditions. Counter ions (Na + or Cl − ) with a concentration equal to 0.9% were included to maintain the neutral state of the systems. To maintain the pressure value at 1 atm, the Berendsen barostat technique [44] was used. The long-range coulomb forces were computed employing the particle-mesh Ewald (PME) method [45,46]. The cut-off radius was set to 8 Å for the non-bonded interactions. The Langevin thermostat method was employed to hold the value of temperature at 300 K [47]. The periodic boundary conditions were also taken into account. The cubic simulation cell was chosen to be larger than the studied proteinligand complexes by 20 Å in every instance. With a multiple time step of 1.25 fs, a regular simulation speed was preserved for intramolecular processes. At an integration step of 2 fs, all intermolecular bonds, including hydrogen bonds, were constrained using the SHAKE algorithm [48]. As a final step, production stages were accomplished over simulation times of 50 ns, 100 ns, 150 ns, and 200 ns. Snapshots of the simulation trajectory were held every 100 ps after an equilibration time of 1-2 ns, determined by the root-mean-square deviations (RMSDs) of the solutes from the initial structure. Simulation steps were executed using a pre-installed macro (md_run.mcr) within the YASARA package. The YASARA program (version 21.11.16) uses the Poisson-Boltzmann approach [49] (called "PBS"). The surfcost parameter, used to calculate the entropic cost of exposing an Å 2 to the solvent, was set to 0.35. The Amber14 force field was used to calculate the binding energy of the inhibitor. To ensure consistency with the empirically determined values, the binding energies derived from PBS were divided by a factor of 20 [50]. The root-mean-square deviation (RMSD), radius of gyration (Rg), root-mean-square fluctuation (RMSF), and solvent-accessible-surface-area (SASA) were all used to analyze the best compounds at the end of the 200 ns MD simulations.

Binding Energy Calculations
Binding free energies of the investigated drugs against Omicron were computed using the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) approach [51]. The below-illustrated equations were used in the process of MM/PBSA binding free energy calculations.
Values of the binding energy of the complex, protein, and ligand are described by ∆G C , ∆G P , and ∆G L , respectively. In addition, ∆G Sol , ∆E MM , and −T∆S stand for the solvation Gibbs energy, gas-phase molecular mechanics change, and conformational entropy, respectively. The term ∆E MM can be determined through summation of the van der Waals and electrostatic interactions. The term ∆G Sol can readily be defined as adding the polar and non-polar solvation values. The entropic contribution is denoted by the term −T∆S.

Drug Likeness Properties
In order to assess the physicochemical parameters of the specified compounds, the online Molinspiration cheminformatics software (http://www.molinspiration.com 30 April 2022) was employed. For each discovered compound, various descriptors were checked, including the octanol/water partition coefficient (milogP), topological polar surface area (TPSA), molecular weight (MWt), number of hydrogen bond donors (nOHNH), number of hydrogen bond acceptors (nON), number of rotatable bonds (Nrotb), and percentage of absorption (%ABS). The equation that was used in computing %ABS is displayed below [52]:

Conclusions
In the current pandemic, with the emergence of Omicron, in silico strategies may be useful in discovering potent inhibitors for this disease. In the present study, molecular docking calculations were performed for seventeen isolated compounds from Euphorbia dendroides. The results showed that six flavonoid compounds are the best anti-Omicron drug candidates, further supporting their efficacy. Combined molecular dynamics simulations (MD) and molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) binding energy evaluations over 200 ns were performed for these six compounds. Two compounds, namely quercetin-3-O-β-D-glucuronopyranoside and quercetin-3-O-glucuronide 6"-O-methyl ester, showed promising binding affinities, with ∆G binding of −41.0 and −32.4 kcal/mol, respectively. As for drug-like properties, both compounds also proved their potential, with a good percentage of absorption (%ABS).