Abstract
The SARS-CoV-2 Omicron variant has spread rapidly and is considered the predominant variant in the world, and its main characteristic is related to evade immunity from natural infection or vaccines, due to its multiple mutations in the spike protein. On the other hand, medicinal plants have been used as alternatives therapies to ameliorate some signs and symptoms in COVID-19, and in our previous work, the cat’s claw (Uncaria tomentosa) stem bark has been studied in vitro and showed antiviral activity on SARS-CoV-2 as well as in silico studies on the 3CLpro protein and as disruptor between the ACE-2 human receptor and the spike protein. The aim in this computational study was to determine the main phytochemical constituents from U. tomentosa stem bark against the SARS-CoV-2 Omicron spike protein based on molecular modeling. A molecular docking was carried out on the isolated phytochemicals in a previous work against the SARS-CoV-2 Omicron spike protein-binding domain (PDB ID: 7T9K). Next, a molecular dynamic study was carried out to monitor the stability during the MD simulations. As results proanthocyanidin-C1 (-10.76 kcal/mol), quinovic acid-type 2 (-9.86 kcal/mol), and proanthocyanidin-B2 (-9.82 kcal/mol) were the constituents with the best binding free energy on the SARS-CoV-2 Omicron spike protein, and the best compound was stable during the dynamic simulation under physiological conditions. It is concluded that the anthocyanidin-based compounds determined in the stem bark ethanol extract could be responsible for the potential antiviral activity on SARS-CoV-2 Omicron variant, and the proanthocyanidin-C1 emerged as a powerful candidate to combat new variants.
1. Introduction
Since November 2021, the World Health Organization (WHO) declared the variant of COVID-19 named Omicron (B.1.1.529) as a variant of concern after its detection for the first time in South Africa [1]. S. S. A. Karim and Q. A. Karim [2] demonstrated that the Omicron variant was characterized by the highest number of spike (S) protein mutations compared to Alpha, Beta, Gamma, and Delta variant of concerns. The fifth variant of SARS-CoV-2 has become dominant in countries Africa, Asia, North America, the United States, and Europe [3]. Indeed, on December 3, the number of confirmed Omicron cases was 219, while on December 9, 2021, 2152 people were affected, which a large number of this variant was described in Denmark, the United Kingdom, and South Africa [4]. The Office for National Statistics (ONS) of the United Kingdom showed that in the week ending January 14, 2022, with omicron as the dominant variant, 1382 deaths were registered in England and Wales [5]. According to Omicron variant, the spike protein (S protein) presented 30 mutations in this protein which were identified such as A67V, H69, V70, T95I, G142D, E156, F157, R158G, G339D, S371L, S373P, S375F, Q493R, G496S, Q498R, L452R, T478K, N501Y, Y505H, T547K, D614G, H655Y, N679K, P681R, N764K, D796Y, N856K, Q954H, N969K,and L981F, meanwhile, the Delta variant ((B.1.617.2) presented mutations at T19R, T95I, G142D, E156,F157, R158G, L452R, T478K, D614G, P681R, P812R, and D950N [6].
On the one hand, the vaccine successes against omicron variant remain unclear compared to their power versus infection with SARS-CoV-2 [3]. Additionally, the emergence of new variants of concern (COV) causes the collapse of the health system especially in developing countries which are unable to ensure the health care offer to their population [7]. Thanks to their effectiveness against viral disease, the plant extracts and their compounds have received increasing attention and might be a hopeful alternative to discover antiviral agents [8–10]. In this context, several research works have attempted to investigate the effect of the plant extracts and natural products against COVID-19. Indeed, in silico study, Alamri and collaborators [11] identified 9 natural compounds including calonysterone, isocodonocarpine, and withanolide A bound strictly to PLpro, while chrysophanol 8-(6-galloylglucoside), luteolin 7-rutinoside, and kaempferol 7-(6″-galloylglucoside) linked strongly with RdRp, and 3,4,5-tri-O-galloylquinic acid, mulberrofuran G, and chrysophanol 8-(6-galloylglucoside) bound efficiently to 3CLpro. In another study, Qamar et al., [12] using in silico approach demonstrated that licoleafol, Calceolarioside B, methyl rosmarinate, myricitrin, and 5,7,3,4-tetrahydroxy-2-(3,3-dimethylallyl) isoflavone from Glycyrrhiza uralensis, Hyptis atrorubens Poit, Myrica cerifera, and Psorothamnus arborescens, respectively, exert an inhibition of SARS-CoV-2 3CLpro effect. These products showed better linking affinity and better docking scores compared to the nelfinavir and prulifloxacin a control drug. Similarly, Yepes-Pérez et al. [13] reported the potential therapeutic activities, in silico, of cadambine, speciophylline, and proanthocyanidin-B2 against SARS-CoV-2 by strong interaction with 3CLpro. Another study showed the inhibition effect of angiotensin-converting enzyme 2 (ACE2) receptor caused by isothymol isolated from Ammoides verticillata plant [14]. A study conducted, in vitro, by Yepes-Pérez and colleagues found an interesting inhibition of 92.7% of SARS-CoV-2 at 25.0 μg/mL and an important reduction the cytopathic effect induced by SARS-CoV-2, after 48 h hydroalcoholic extract of Uncaria tomentosa treatment on Vero E6 cells [15]. In recent study, Nair et al. [16] investigated, in vivo, the antiviral activity of Artemisia annua L. hot-water extracts on five infectious variants of SARS-CoV-2 including alpha (B.1.1.7), beta (B.1.351), gamma (P.1), delta (B.1.617.2), and kappa (B.1.617.1) using Vero E6 cells. As results, the extracts of A. annua exhibit an important activity versus SARS-CoV-2 and its five variants.
Uncaria tomentosa, a member of Rubiaceae family, commonly called cat’s claw or uña de gato is distributed in the Amazon region of Peru and other South American countries such as Colombia, Ecuador, and Brazil [17]. U. tomentosa commonly has been used by patients to treat the gastrointestinal disease, inflammatory diseases, and viral infections [18]. Previous papers reported that U. tomentosa contains various secondary metabolites including triterpenes, oxindole alkaloids, vegetal steroids, tannin, phenolic compounds, and flavonoids [19, 20]. Concerning the biological proprieties, pharmacological investigations showed that plant extracts exhibit multiple biological properties as antiviral on RNA-type virus like dengue [21], herpes [22], and SARS-CoV-2 [15]. Molecular docking protocols are a valuable computational strategy used in modern drug design which studies the ability of one or a few compounds to bind a given molecular target. In this regard, protein-ligand dockings of the main components of ethanolic extracts of cat’s claw were performed against the three-dimensional structure of the SARS-CoV-2 Omicron spike glycoprotein (PDB ID: 7T9K) [23] (Figure 1).
2. Materials and Methods
2.1. Ligand Dataset Preparation of U. tomentosa Stem Bark
Ligands used in this computational study are components determined of the cat’s claw ethanol extract by HPLC-MS in the antiviral activity in vitro against SARS-CoV-2. These compounds were classified as oxindole alkaloids: isomitraphylline, isopteropodine, isorynchophylline, mitraphylline, pteropodine, rhynchophylline, speciophylline, and uncarine F; indole alkaloid: 3-dihydrocadambine; quinovic acids: quinovic acid-2; and proanthocyanins: chlorogenic acid, epicatechin, proanthocyanidin-B2, proanthocyanidin-B4, and proanthocyanidin-C1.
2.2. Generation of 3D-Pharmacophore Model
SARS-CoV-2 Omicron spike is one of the most essential and promising targets; it is the main part that SARS-CoV-2 virus used to penetrate the host cells. Spike is consisted of two main parts S1 and S2; S1 part contain receptor-binding domain (RBD); it can recognize and bind with ACE-2 that facilitate the attack of virus against targeted host cells Figure [2]. Azithromycin shows a good inhibitory activity to RBD that can restrict the binding of SARS-CoV-2 spike with host ACE-2 [24]. The structure-based pharmacophore model was generated by BIOVIA Discovery Studio 2019 client software; the pharmacophore queries was generated depending on the amino acid present in RBD; the best-fitted values and relative fit values were calculated in Table 1. The main objective of the pharmacophore structure-based drug design is to generate a 3D-pharmacophore model based on the known pocket. In this study, the pocket of RBD-SARS-CoV-2 was used to generate the queries (features) based on the amino acids present in target pocket. These features (Figures 2, 3, and 4) are essential for the activity. Consequently, any ligand will match or fit with these features is expected to be active.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
2.3. Molecular Docking Studies
After doing the structure-based pharmacophore study for fifteen compounds isolated from Uncaria tomentosa (cat’s Claw) stem bark, the molecular docking was done inside the pocket of RBD-SARS-CoV-2 Omicron spike protein by using BIOVIA Discovery Studio 2019 client. The SARS-CoV-2 Omicron spike protein crystal protein (PDB code: 7T9K) was downloaded from https://www.rcsb.org, acceded on January 25, 2022. At first water molecules have been removed from the complex. Then, crystallographic disorders and unfilled valence atoms were corrected using protein report and utility and clean protein options. Protein energy was minimized by applying CHARMM force fields. The rigid of binding site was the structure of protein that was obtained by applying fixed atom constraint. The protein essential amino acids were defined and prepared for docking process. 2D structures of tested compounds were drawn using ChemBioDraw Ultra17.0 and saved in MDL-SD file format; from BIOVIA Discovery Studio 2019 client, the saved file was opened, 3D structures were protonated, and energy was minimized by applying 0.05 RMSD kcal/mol CHARMM force field. Then, the minimized structures were prepared for docking using prepared ligand protocol. The molecular docking process was carried out using the CDOCKER protocol. The receptor was held rigid, while the ligands were allowed to be flexible during the refinement in which each molecule was allowed to produce ten different interaction poses with the protein. Then docking scores (-CDOCKER interaction energy) of the best-fitted poses with the active site at RBD-SARS-CoV-2 Omicron spike protein were recorded, and a 3D view was generated [25]. We use all these processes to predict the proposed binding mode, affinity, preferred orientation of each docking pose, and binding free energy (∆G) of the tested compounds with RBD-SARS-CoV-2 Omicron spike protein Table 2.
2.3.1. Validation of Molecular Docking
The molecular docking was initially validated by redocking of the cocrystallized ligand into the active site of the respective receptor with the calculation of root mean square deviation (RMSD) for reliability and reproducibility of the proposed docking algorithm. In cocrystallized ligand not integrated with SARS-CoV-2 Omicron spike, we used Azithromycin as a guide and as coligand. Azithromycin was redocked on RBD-SARS-CoV-2 Omicron spike protein (PDB ID: 7T9K) and the RMSD 1.46 blew 2.00 A° indicating a validated method.
2.4. Molecular Dynamics Simulation (MD)
Molecular dynamics simulation of the protein-ligand complexes was carried on the docked complexes for proanthocyanidin-B2 (the most active in the molecular docking analysis) with the SARS-CoV-2 Omicron spike protein (PDB ID: 7T9K) using GROMACS [26] 2021.1 version and Linux 5.4 package. The GROMOS96 54a7 force field was selected as the force field for proteins, and the ligand topologies were generated from the PRODRG [27] server. All the complexes were solvated using simple point charge (SPC) water molecules in a rectangular box. To make the simulation system electrically neutral, the required number of Na+ and Cl− ions were added, while 0.15 mol/L salt concentrations were set in all the systems. Using the steepest descent method, all the solvated systems were subjected to energy minimization for 5000 steps. Afterwards, NVT (constant number of particles, volume, and temperature) series, NPT (constant number of particles, pressure, and temperature) series, and the production run were conducted in the MD simulation. The NVT and the NPT series were conducted at a 300 K temperature and 1 atm pressure for the duration of 300 ps. V-rescale thermostat and Parrinello-Rahman barostat were selected for the performed simulation. Finally, the production run was performed at 300 K for a duration of 100 ns (nanoseconds). Thereafter, a comparative analysis was performed measuring root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and hydrogen bonds to analyze their stability. The Xmgrace program was used to represent the analyses in the form of plots.
2.5. Molecular Similarity
The molecular similarity is a computational ligand-based calculation that is used to expect the similarity between two molecules. The principle that this study based on that the structurally similar ligands are predicted to have related properties to different molecules. Molecular similarity calculation targets the molecular properties of the examined ligands, such as LogP and molecular weight. Additionally, the molecular similarity mostly depends on the measures of distances in space in both structural and physicochemical descriptors.
3. Results and Discussion
3.1. Selecting of Ligands Based on Phytochemical Constituents of U. tomentosa Ethanol Extract
The ligands selected were established based on our previous study, which 15 metabolites were determined by HPLC-MS [15] in hydroalcoholic extract of U. tomentosa (Figure 2). They were identified as oxindole alkaloids, indole alkaloids, quinovic acid, and proanthocyanidins. The total phytoconstituents were classified as spirooxindole alkaloids, indole glycosides alkaloids, quinovic acid, and proanthocyanidins. Other investigations have confirmed our findings based on the phytochemical profile carried out on different extracts from stem bark, which is traditionally consumed in different forms such as infusions, teas, decoction, macerated, and lately dosage forms like capsules and tablets which can be sold in pharmacies. Although in national hospitals from Peru, bags containing a powder of U. tomentosa (uña de gato) are used for rheumatic pain as alternative treatment. Several studies attributed to its alkaloids as the responsible of the biological activity like antiviral. However, other studies revealed the total extract synergisms of the pharmacological effect, acting by different mechanisms. Furthermore, the quality control of U. tomentosa is based on the determination of its eight alkaloids as in the American Pharmacopoeia (USP) and corresponds to its fingerprint [28].
3.2. Generation of 3D-Pharmacophore Model
3.2.1. Activity Prediction
The test set of the isolated compounds and Azithromycin were examined against the generated 3D-pharmacophore model. In this process, fit values and relative fit values (that quantitatively represent the existence of that essential features) were calculated. All ligands were proved to have the main essential features of inhibitor (Table 1). Interestingly, some ligands showed high fit values and good matching with 3D-pharmacophore queries.
3.3. Docking Studies
Molecular docking protocol is a valuable computational strategy used in modern drug design which study the ability of one or a few compounds to bind a given molecular target. In this regard, protein-ligand dockings of the main components of ethanolic extracts of cat’s claw were performed against the three-dimensional structure Omicron spike glycoprotein-RBD of SARS-CoV-2 (PDB ID: 7T9K) [23].
The binding mode of Azithromycin exhibited an energy binding of -9.32 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which formed nine Pi-Alkyl interactions with Val976, Met740, Leu966, Pro589, and Lys856; moreover, four hydrogen bonds was created with Thr573, Thr549, and Cys743 with a distance of 1.78, 1.82, 1.84, and 1.87 Å, respectively, additionally the cationic NH2 center binding with Asp745 ionic interaction (Figure 5).
The binding mode of 3-dihydrocadambine exhibited an energy binding of -8.58 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which interacted with Lys547, Cys590, Gly744, Asp745, Thr549, and Lys856 by eleven hydrogen bonds with a distance of 2.72, 2.77, 2.60, 2.50, 2.59, 1.86, 2.16, 2.60, 2.51, 2.29, and 2.68°A, and additionally formed two Pi-Alkyl interactions with Val976 and Lys856 (Figure 6).
The binding mode of chlorogenic acid exhibited an energy binding of -6.42 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which creates one Pi-alkyl interaction with Leu966 and additionally interacts with Lys856, Tyr741, Met740, Lys547, Gly744, Asp745, Thr549, Leu977, and Asn978 by eleven hydrogen bonds with a distance range of 1.91 and 2.90°A (Figure S1). The binding mode of isomitraphylline exhibited an energy binding of -7.75 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which created three Pi-Alkyl interactions with Val976, Lys856, and Pro589; additionally interacted with Leu966, Gly744, Asn978, Asp745, and Thr549 by seven hydrogen bonds; and moreover, interacted with Asp745 by ionic interaction with cationic amino group. (Figure 7).
The binding mode of isopteropodine exhibited an energy binding of -7.66 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. The hydrophilic centers create four conventional hydrogen bonds and seven nonconventional hydrogen bonds with Gly548, Thr549, Asn540, Asp745, Glu748, and Asn978 and additionally interact with Leu977 and Val976 by two Pi-alkyl interactions and form ionic interactions with Glu748 and Asp745 (Figure 8).
The binding mode of isorhynchophylline exhibited an energy binding of -7.35 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which creates five Pi-Alkyl and Pi-Sigma interactions with Ile587, Pro589, Phe541, and Val976. Besides that, it was interacted with Thr572, Thr573, and Lys 856 by four hydrogen bonds with a distance of 2.05 to 2.68°A (Figure 9).
The binding mode of mitraphylline exhibited an energy binding of -7.67 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which creates two Pi-Alkyl interactions with Pro589 and Leu977 and, moreover, interacts with Gly744, Asp745, Thr549, and Lys547 by five hydrogen bonds with a distance range of 2.09 to 2.68°A (Figure 10).
The binding mode of proanthocyanidin-B2 exhibited an energy binding of -9.82 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. This complex was stabilized by creating four strong Pi-Alkyl, Pi-amide, and Pi-sigma interactions with Lys547, Gly548, Pro589 and Thr572 and, moreover, interacted with Arg567, Thr573, Gly744, Asn978, Asp745, Thr549, Ile587, and Thr572 by nine hydrogen bonds with a distance range of 2.00 to 3.07°A (Figure 11A), while proanthocyanidin-B4 exhibited an energy binding of -8.95 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. It creates three Pi-Alkyl and Pi-amide interactions with Cys743, Met740, and Val976 and additionally interacts with Thr573, Cys590, Thr549, Pro589, Asp745, Gly744, Ser975, and Lys547 by eight hydrogen bonds with a distance range of 1.8 to 3.02°A (Figure 11(b)).
(a)
(b)
The binding mode of proanthocyanidin-C1 exhibited an energy binding of -10.76 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. The complex was stabilized by creating four Pi-Alkyl and Pi-Sigma interactions with Lys856, Pro589, and Val976 and, moreover, interacted with Gly744, Pro589, Ser967, Thr549, Asp745, Tr572, and Thr573 by nine hydrogen bonds with a distance range of 2.11 to 2.99°A.(Figure 12).
The binding mode of quinovic acid exhibited an energy binding of -9.86 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which creates three Pi-Alkyl interactions with Lys547 and Val976 and moreover interacts with Asp745, Gln321, Cys743, Arg1000, Tyr741, and Thr573 by seven hydrogen bonds with a distance range of 1.68 to 2.52 Å (Figure 13).
The binding mode of rhynchophylline exhibited an energy binding of -7.49 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. It creates four Pi-Alkyl interactions with Leu966, Met740, Phe592, and Pro589 and, additionally, interacts with Gly458 and Thr549 by two hydrogen bonds with a distance of 2.57 and 2.58 Å (Figure 14). The binding mode of pteropodine exhibited an energy binding of -6.72 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which creates four Pi-Alkyl interactions with Leu546, Phe541, Val976, and Leu966 and moreover interacts with Gly548, Val976, and Ser975 by three hydrogen bonds (Figure S2).
The binding mode of speciophylline exhibited an energy binding of -6.45 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein, which creates three Pi-Alkyl interactions with Leu966, Leu977, and Met740 and moreover interacts with Arg1000 and Thr549 by two hydrogen bonds (Figure S3A), while uncarine F exhibited an energy binding of -7.11 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. It forms two Pi-Alkyl interactions with Pro589 and Leu977 and moreover interacts with Asp745, Gly744, Thr549, and Lys547 by four hydrogen with a distance range of 2.09 to 2.68 Å (Figure S3B). The binding mode of epicatechin exhibited an energy binding of -6.89 kcal/mol-1 against RBD-SARS-CoV-2 Omicron spike protein. It creates one Pi-Alkyl interaction with Pro589 and additionally interacts with Ile578, Lys856, and Asp745 by three hydrogen bonds (Figure S4).
These facts are in good agreement with our previous report which computational and in vitro experiments confirmed that the hydroalcoholic fraction from cat’s claw inhibited the release of SARS-CoV-2 infectious particles and reduced the cytopathic effect caused by the virus [15]. Therefore, virtual screening suggests that some components in hydroalcoholic extract of U. tomentosa stem bark could positively modulates the recently emerged SARS-CoV-2 Omicron infection. However, experimental assays are needed to support these computational predictions. In addition, as can be seen in a previous report [28], calculated pharmacokinetic indices for the most qualified cat’s claw components showed the druggability of the selected components and their potential as likely orally active antiviral. Calculated parameters revealed that most of components exhibited favorable characteristics as the drug like; hence, U. tomentosa may constitute itself a promissory option to fight against COVID-19 infection.
3.4. Molecular Dynamics Simulation (MD)
Molecular dynamics and simulation (MD) studies were carried out in order to determine the stability and convergence of 7T9K/proanthocyanidin-B2, 7T9K/proanthocyanidin-C1, and 7T9K/quinovic acid complexes. According to 7T9K/Pproanthocyanidin-B2, the dynamic movements of atoms and conformational variations of backbone atoms of the protein-ligand complex were calculated by RMSD to detect their stability upon apo and ligand bonded state. It is observed that the protein and ligand exhibit very lower RMSD with no major fluctuations indicating their greater stability. The complex was slightly stable until 25 ns~ and showed more stability later. The flexibility of each residue was calculated in terms of RMSF to get better insight on the region of proteins that are being fluctuated during the simulation. It can be understood that the binding of ligand does not make the protein in any areas. The compactness of the complex was represented by the radius of gyration (Rg). The lower degree of fluctuation throughout the simulation period indicates the greater compactness of a system. The Rg of the complex was found to be slightly lower than the starting period. Interaction between protein-ligand complexes and solvents was measured by solvent accessible surface area (SASA) over the simulation period. So, SASA of the complex was calculated to analyze the extent of the conformational changes occurred during the interaction. Interestingly, the protein featured a slight reduction of the surface area showing relatively lower SASA value than the starting period. Hydrogen bonding between a protein-ligand complex is essential to stabilize the structure. It was observed that the highest number of conformations of the protein formed up to six hydrogen bonds with the ligand (Figures 15, 16, 17, 18, 19, and 20). Additionally, we calculated the binding free energy of the last 20 ns of MD production run of the protein-ligand complex with an interval of 100 ps from MD trajectories using MM/PBSA method. We also utilized the MmPbSaStat.py script that calculated the average free binding energy and its standard deviation/error from the output files that were obtained from g_mmpbsa. The ligand showed less binding free energy of -121 KJ/mol with the protein. Further, we identified the contribution of each residue of the protein in terms of binding free energy to the interaction with the ligand. By decomposing the total binding free energy of the system into per residue contribution energy, the contribution of each residue was calculated. This gave us an insight into the “crucial” residues that contributes favorably to the binding of this molecule to the protein. It was found that THR-549, ARG-567, THR-573, ILE587, GLY-744, ASP745, and ASN-978 residues of the protein contributed higher than -5 KJ/mol binding energy and thereby are hotspot residues in binding with the ligand.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
According to 7T9K/proanthocyanidin-C1 complex, it is observed that the protein and ligand exhibit very lower RMSD with no major fluctuations indicating their greater stability. The complex was unstable until 10 ns~, and it becomes stable until 40 ns~ with minor fluctuation and showed stability later. The flexibility of each residue was calculated in terms of RMSF to get better insight on the region of proteins that are being fluctuated during the simulation. It can be understood that the binding of ligand makes the protein slightly flexible in 1700-1800 residue areas. The compactness of the complex was represented by the radius of gyration (Rg). The lower degree of fluctuation throughout the simulation period indicates the greater compactness of a system. The Rg of the complex was found to be slightly lower than the starting period. The interaction between protein-ligand complexes and solvents was measured by solvent accessible surface area (SASA) over the simulation period. So, SASA of the complex was calculated to analyze the extent of the conformational changes occurred during the interaction. Interestingly, the protein featured a slight reduction of the surface area showing relatively lower SASA value than the starting period. Hydrogen bonding between a protein-ligand complex is essential to stabilize the structure. It was observed that the highest number of conformations of the protein formed up to two hydrogen bonds with the ligand. Moreover, we calculated the binding free energy of the last 20 ns of MD production run of the protein-ligand complex with an interval of 100 ps from MD trajectories using MM/PBSA method. We also utilized the MmPbSaStat.py script that calculated the average free binding energy and its standard deviation/error from the output files that were obtained from g_mmpbsa. The ligand showed less binding free energy of -96 KJ/mol with the protein. Further, we identified the contribution of each residue of the protein in terms of binding free energy to the interaction with the ligand. By decomposing the total binding free energy of the system into per residue contribution energy, the contribution of each residue was calculated. This gave us an insight into the “crucial” residues that contributes favorably to the binding of this molecule to the protein. It was found that LYS-547, THR573, TYR-741, ASP745, GLU-784, and ARG-1000 residues of the protein contributed higher than -8 KJ/mol binding energy and thereby are hotspot residues in binding with the ligand.
According to 7T9K/quinovic acid complex, the dynamic movements of atoms and conformational variations of backbone atoms of the protein-ligand complex were calculated by RMSD to detect their stability upon apo and ligand bonded state. It is observed that the protein, ligand, and complex exhibit very lower RMSD with no major fluctuations indicating their greater stability. The complex was unstable until 17 ns~ and showed stability later. The flexibility of each residue was calculated in terms of RMSF to get better insight on the region of proteins that are being fluctuated during the simulation. The compactness of the complex was represented by the radius of gyration (Rg). The lower degree of fluctuation throughout the simulation period indicates the greater compactness of a system. The Rg of the complex was found to be slightly lower than the starting period. Interaction between protein-ligand complexes and solvents was measured by solvent accessible surface area (SASA) over the simulation period. So, SASA of the complex was calculated to analyze the extent of the conformational changes occurred during the interaction. Interestingly, the protein featured a slight reduction of the surface area showing relatively lower SASA value than the starting period. Hydrogen bonding between a protein-ligand complex is essential to stabilize the structure. It was observed that the highest number of conformations of the protein formed up to five hydrogen bonds with the ligand.
Additionally, we calculated the binding free energy of the last 20 ns of MD production run of the protein-ligand complex with an interval of 100 ps from MD trajectories using MM/PBSA method. We also utilized the MmPbSaStat.py script that calculated the average free binding energy and its standard deviation/error from the output files that were obtained from g_mmpbsa. The ligand showed less binding free energy of -120 KJ/mol with the protein. Further, we identified the contribution of each residue of the protein in terms of binding free energy to the interaction with the ligand. By decomposing the total binding free energy of the system into per residue contribution energy, the contribution of each residue was calculated. This gave us an insight into the “crucial” residues that contributes favorably to the binding of this molecule to the protein. It was found that THR-549, THR-572, GLY-744, and ASP745 residues of the protein contributed higher than -7 KJ/mol binding energy and thereby are hotspot residues in binding with the ligand.
Finally, the MD simulation results suggest that (1) the ligands were stable at the binding pocket in the 100-ns MD simulations, (2) ligand does not leave the Omicron receptor-binding domain (RBD) while 100 ns MD run, and (3) 100 ns molecular dynamics simulation suggest the rationality and validity of our docking studies. Accordingly, some chemical components of hydroalcoholic extract of U. tomentosa, particularly, proanthocyanidin-B2, proanthocyanidin-C1, and quinovic acid could suppress Omicron viral spike function, which in turn suggest that U. tomentosa may interfere with internalization of the Omicron variant in the human host.
Altogether, computational modeling studies showed strong evidence that several chemical components in hydroalcoholic extract of U. tomentosa, especially proanthocyanidin-B2, proanthocyanidin-C1, and quinovic acid are able to bind to SARS-CoV-2 Omicron spike protein via key interactions with those critical amino acids which mediates attachment of the virus to the host cell surface in stable complexes. In turn, we firmly believe that these results could open new avenues for treating COVID-19 infection associated with the Omicron variant. However, further studies are needed aiming to validate, support, and extend computational findings, as well as determine the potential benefits of hydroalcoholic extract of U. tomentosa against the SARS-CoV-2 Omicron variant.
3.5. Molecular Similarity
In this research, a molecular similarity study has been preceded for fifteen ligands that showed antiviral activities against Azithromycin targeted on RBD-SARS-CoV-2 Omicron spike protein, by using Discovery studio software. The used molecular properties listed in Table 3 [29] include the number of rotatable bonds, number of cyclic rings, number of aromatic rings, number of hydrogen bond donors (HBD), number of hydrogen bond acceptors (HBA), partition coefficient (ALog p), molecular weight (M. Wt), and molecular fractional polar surface area (MFPSA) (Figure 21) which shows that rhynchophylline, isorhynchophylline, proanthocyanidin-B4, proanthocyanidin-B2, and quinovic acid (red balls) are similar to Azithromycin (green ball).
4. Conclusions
U. tomentosa have demonstrated antiviral activity against the SARS-CoV in vitro and also in silico on potential targets such as the RBD in spike protein, the complex SARS-CoV-2/ACE-2 junction, and the 3CLpro protein., Based on the phytochemicals constituents determined in our previous work, 15 compounds were docked against SARS-CoV-2 Omicron spike protein-binding domain (RBD) (PDB ID: 7T9K) resulting proanthocyanidin-C1 (-10.76 Kcal/mol), quinovic acid-type 2 (-9.86 kcal/mol), and proanthocyanidin-B2 (-9.82 kcal/mol) like the constituents with best binding energy on the SARS-CoV-2 omicron spike protein and proanthocyanidin-C1 which was stable at the binding pocket in the 100-ns of molecular dynamic simulations under physiological parameters. Although in vitro analysis on omicron variant is needed, proanthocyanidin-C1 might be considered as a candidate molecule against new variants or the total extract of U. tomentosa, which all constituents would synergism the antiviral effect.
Data Availability
The data used to support the findings of this study are included within the supplementary information file.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
Conceptualization was contributed by O.H.-C., A.M.S., and A.F.Y.-P.; methodology was contributed by N.H.A.; software was contributed by S.A. and H.C.; validation was performed by A.M.S., A.F.Y.-P., E.L.-G., R.D.H.-Q., and J.S.A.-G.; formal analysis was performed by A.M.S. and G.E.-S.B.; investigation was performed by O.H.-C. and J.F.K.-C.; writing—original draft preparation was initiated by O.H.-C., H.C., and J.B.P.-O.; writing—review and editing was initiated by A.M.S. and T.B.; project administration was performed by S.A.; funding acquisition was performed by N.H.A. All authors have read and agreed to the published version of the manuscript.
Acknowledgments
This work was funded by Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2022R62), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia, and King Saud University Researchers Supporting Project number (RSP-2021/26), Riyadh, Saudi Arabia.
Supplementary Materials
Figure S1: chlorogenic acid docked in RBD-SARS-CoV-2 Omicron spike protein, hydrogen bonds (green), and the pi interactions are represented in purple lines. Figure S2: pteropodine docked in RBD-SARS-CoV-2 Omicron spike protein, hydrogen bonds (green), and the pi interactions are represented in purple lines. Figure S3: (A) speciophylline docked in RBD-SARS-CoV-2 Omicron spike protein and (B) uncarine F docked in RBD-SARS-CoV-2 Omicron spike protein, hydrogen bonds (green), and the pi interactions are represented in purple lines. Figure S4: epicatechin docked in RBD-SARS-CoV-2 Omicron spike protein, hydrogen bonds (green), and the pi interactions are represented in purple lines. (Supplementary Materials)