Structure-based virtual screening and molecular dynamics of phytochemicals derived from Saudi medicinal plants to identify potential COVID-19 therapeutics.

Coronavirus disease 2019 (COVID-19) has affected almost every country in the world by causing a global pandemic with a high mortality rate. Lack of an effective vaccine and/or antiviral drugs against SARS-CoV-2, the causative agent, has severely hampered the response to this novel coronavirus. Natural products have long been used in traditional medicines to treat various diseases, and purified phytochemicals from medicinal plants provide a valuable scaffold for the discovery of new drug leads. In the present study, we performed a computational screening of an in-house database composed of ~1000 phytochemicals derived from traditional Saudi medicinal plants with recognised antiviral activity. Structure-based virtual screening was carried out against three druggable SARS-CoV-2 targets, viral RNA-dependent RNA polymerase (RdRp), 3-chymotrypsin-like cysteine protease (3CLpro) and papain like protease (PLpro) to identify putative inhibitors that could facilitate the development of potential anti-COVID-19 drug candidates. Computational analyses identified three compounds inhibiting each target, with binding affinity scores ranging from -9.9 to -6.5 kcal/mol. Among these, luteolin 7-rutinoside, chrysophanol 8-(6-galloylglucoside) and kaempferol 7-(6″-galloylglucoside) bound efficiently to RdRp, while chrysophanol 8-(6-galloylglucoside), 3,4,5-tri-O-galloylquinic acid and mulberrofuran G interacted strongly with 3CLpro, and withanolide A, isocodonocarpine and calonysterone bound tightly to PLpro. These potential drug candidates will be subjected to further in vitro and in vivo studies and may assist the development of effective anti-COVID-19 drugs.

Introduction

The global threat of the 2019 novel coronavirus disease (now officially named COVID-19 and the causative pathogen is SARS-CoV-2) is rapidly escalating with unprecedented international health and economic burden in the recent history. The whole-genome of SARS-CoV-2 had been sequenced and revealed that SARS-CoV-2 pathogen is the fifth strain of β-coronaviruses, which include OC43, HKU1, SARS-CoV-1 and Middle East respiratory syndrome coronavirus (MERS-CoV) (Wang et al., 2020). After the entry into the host cell through strong binding of β-coronaviruses’ protein spikes with the human angiotensin-converting enzyme 2 (ACE2) receptor (Xu et al., 2020), β-coronaviruses generate large polyproteins (PP1a and PP1ab) upon genome translation of open reading frame (ORF) 1a and ORF1ab by the host cell machinery. These polyproteins, also known as replicase polyproteins, are proteolytically cleaved by essential cysteine proteases encoded by the virus, explicitly 3-chymotrypsin-like protease (3CLpro; sometimes called main protease (Mpro)) and papain-like proteases (PLpro), to release 16 non-structural proteins (nsps) (Thiel et al., 2003, Ziebuhr, 2004). Among those nsps is the RNA-dependent RNA polymerase (RdRp or sometimes referred to as nsp12). RdRp catalyses the synthesis of a complementary RNA strand using the viral RNA and hence, it plays an essential role in directing the replication and transcription of SARS-CoV-2 genome (Chen et al., 2020, Lung et al., 2020, Zumla et al., 2016). The proteolysis of PP1a and PP1ab by 3CLpro occurs at 11 distinct sites and generates various nsps that are important for the viral replication (Anand et al., 2003). Therefore, 3CLpro  is crucial for the virus particle replication and represent a valid target for the identification of coronavirus inhibitors (Needle et al., 2015). PLpro is another crucial cysteine protease, an enzyme that cleaves N-terminus of the replicase polyprotein to release several nsps, among them the nsp3, in which PLpro is encoded and is implicated not only in the viral replication but also in suppressing the host innate immune response (Chen et al., 2020, Lei et al., 2018, Yuan et al., 2015). In addition, PLpro possesses a nucleic acid-binding domain (NAB) with a nucleic acid chaperon function and is essential in the virus replication correction. Furthermore, PLpro has posttranslational modification activities, deubiquitination and deISGylating, on the host proteins, such as interferon regulatory factor 3(IRF3), which consequently lead to the removal of ubiquitin and human interferon-stimulated gene 15 (ISG15), the main signalling elements of the antiviral innate immune response (Báez-Santos et al., 2015, Chen et al., 2020, Lei et al., 2018, Liu et al., 2020). The indispensable role of the RdRp, 3CLpro and PLpro in orchestrating the RNA genome replication/transcription cycles of SARS-CoV-2 make them potential targets for designing antiviral candidates against COVID-19 (Báez-Santos et al., 2015, Ramajayam et al., 2011, Ren et al., 2013).

For ages, phytochemicals have been found to be a fruitful source of molecules with diverse therapeutic potentials and still are considered a valuable resource for the discovery of novel drug leads in other word, herbal medicines and purified phytochemicals can be used to develop more efficient drugs based on the structure of natural compounds (Mani et al., 2020). Given the time of its emergence, few studies related to the development of naturally-derived inhibitors of three main druggable targets of SARS-CoV-2 (RdRp, 3CLpro and PLpro) have already been reported using computer modelling for screening purposes. For instance, in silico or biological screening of a series of biologically active natural compounds have been demonstrated to directly inhibit these important proteins in pervious HCoVs such as SARS-CoV, MERS-CoV (Park et al., 2012, Shen et al., 2019) and the novel SARS-CoV-2 (Aanouz et al., 2020, Ul Qamar et al., 2020, Zhang et al., 2020). (Lung et al., 2020) virtually screened 83 compounds found in Chinese traditional medicines with potential activity against SARS-CoV-2 and identified theaflavin as a potential inhibitor of the RNA-dependent RNA polymerase. Likewise, (Zhang et al., 2020) screened 115 compounds that are commonly used in treating viral respiratory infection in China and highlighted 13 compounds for further studies to prove their activities against SARS-CoV-2. A virtual screening and molecular modelling study of marine natural products have also been investigated and proposed seventeen molecules with potent inhibitory action against 3CLpro (Gentile et al., 2020). In the present study, we conducted a computational study on a library of phytochemicals from Saudi medicinal plants with potential antiviral activities (Aati et al., 2019, Arbab et al., 2017) to identify potent inhibitors of RdRp, 3CLpro and PLpro, which eventually can be employed as novel lead molecules to accelerate the pace of designing anti-COVID-19 candidates.

Materials and methods

Preparation of protein structures

The X-ray crystal structures of three SARS-CoV-2 proteins including the RNA-dependent RNA polymerase (RdRp) with 2.50 Å resolution (accession number 7BV2) (Yin et al., 2020), 3-chymotrypsin-like cysteine protease (3CLpro) with 2.1 Å resolution (accession number 6 W63) (Mesecar, 2020) and papain like protease (PLpro) with 2.7 Å resolution (accession number 6W9C) (Osipiuk et al., 2020) were obtained from the Protein Data Bank. Discovery studio program (Studio, 2008) was used to prepare the proteins for virtual screening and molecular docking by removing co-crystalized ligands and water molecules and to prepare the proteins in monomer form in PDB formats.

Phytochemicals database

A library containing ~1000 plant-derived compounds obtained from 60 traditional Saudi medicinal plants with antiviral activity were designed. The phytochemicals were collected by comprehensive literature study. The 3D structures of compounds were obtained either from PubChem database (available online: https://pubchem.ncbi.nlm.nih.gov/) or were drown manually using ChemDraw Ultra 7.0 and saved in one SDF file.

Structure-based virtual screening

The virtual screening of phytochemical database against the structure of target proteins was performed individually using Autodock vina in PyRx 8.0 virtual screening tool (Dallakyan and Olson, 2015). At first, compounds were imported into OpenBabel program (Babel, 2011) implemented in PyRx for energy minimization using MMFF94 force field and for converting compounds into Autodock PDBQT format. The 3D grid box parameters were set to cover the active site cavity within each protein. The top three compounds with highest binding energy scores against each viral protein were selected as hit compounds for further assessment.

Molecular docking protocol

The molecular docking of candidate phytochemicals against their viral protein’s targets were carried out using Autodock vina 1.1.2 program (Trott and Olson, 2010). Autodock Tools 1.5.6 program (Huey and Morris, 2008) was employed to prepare the PDB structures of proteins for docking by adding polar hydrogen and to convert proteins and ligands PDB files into Autodock PDBQT format. The grid box parameters for RdRp were set to size 24 Å × 22 Å × 28 Å (x, y and z) and centre 113.78 Å × 118.474 Å × 132.008 Å (x, y and z). For 3CLpro, the grid box parameters were set to size 18 Å × 18 Å × 18 Å (x, y and z) and centre −20.318 Å × 19.141 Å × −28.343 Å (x, y and z). As for PLpro, the grid box parameters were set to size 21 Å × 22 Å × 22 Å (x, y and z) and centre −27.423 Å × 30.003 Å × 27.531 Å (x, y and z). The 2D ligand-protein interaction diagrams were generated by PoseView tool (available online: https://proteins.plus/). PyMOL (DeLano, 2002) was used to visualize and analyse the docking results.

Molecular dynamic simulation

All atom molecular dynamic (MD) simulation of ligand–protein complexes were performed with GROMACS 2018 software OPLS-AA/L force field. The high-resolution X-ray structures of proteins (PDB IDs: 6M71, 6W63 and 6W9C for RdRp, 3CLpro and PLpro, respectively) were used as initial structures for the MD simulation. The protein structures were further optimized for simulation by DockPerp tool in UCSF Chimera (Pettersen et al., 2004). The top-predicted docking pose of each ligand with the highest binding score was used as a starting point for the MD simulation. The ligands were parametrized using SwissParam webserver (Zoete et al., 2011). All systems were simulated for 20 ns using previously reported protocol (Alamri, 2020). Basic MD simulation parameters including real-mean-square-deviations (RMSD), real-mean-square-fluctuation (RMSF), radius of gyration (Needle et al., 2015), number of hydrogen bonds (H-bonds) were analysed for each complex.

Results and discussion

Construction of potential phytochemical database

Initially, an in-house phytochemical database composed of ~1000 bioactive compounds derived from 60 traditional Saudi medicinal plants was designed. The extracts of the selected plants were reported to possess antiviral activity against a wide-range of viruses such as HIV, HCV and HCB at various concentrations (Aati et al., 2019, Arbab et al., 2017). The main constituents of each plant were identified by comprehensive literature study. This approach was applied to build a focused phytochemicals library with recognised antiviral activity, so it can be further screened in search for potential anti-COVID19 agents. The tested library composed of phytochemicals that represent major classes of natural products as shown in (Fig. 1 ). The natural occurrence and examples of the top five chemical classes are presented in Table S1.

Fig. 1

The major classes of natural products in the tested library.

Structure-based virtual screening, molecular docking and MD simulation

A high-throughput virtual screening of a diverse plant-derived compounds from Saudi medicinal plants, followed by molecular docking and all atom 20 ns MD simulation studies were performed. The screening was conducted against three well-established drug discovery targets including SARS-CoV-2 RdRp, 3CLpro and PLpro enzymes.

RNA-dependent RNA polymerase (RdRp)

First, we screened the designed database against SARS-CoV-2 RdRp. The RdRp is a conserved enzyme in coronaviruses and essential for their replication and transcription of their genomes (Chen et al., 2020, Lung et al., 2020, Zumla et al., 2016). RdRp in its active form tends to form a heterotrimeric complex with nsp7 and nsp8 for maximum activity (Gao et al., 2020). The virtual screening of in-house phytochemical database was performed against the recently solved X-ray structure of the nsp12 protein. The analysis resulted in the identification of three compounds, namely chrysophanol 8-(6-galloylglucoside), luteolin 7-rutinoside and kaempferol 7-(6″-galloylglucoside), which have been reported to be isolated form Rumex dentatus (Mishra et al., 2018), Marrubium vulgare (Neamah et al., 2009) and Acacia species (El-Toumy et al., 2018), respectively. The natural source and docking scores of these compounds are presented in (Table 1 ).

Table 1

Potential inhibitors of RdRp from in-house Saudi medicinal database.

Name/ID	Chemical structure	Natural source	Docking score (kcal/mol)
Remdesivir monophosphate		–	−7.0
Chrysophanol 8-(6-galloylglucoside)(PubChem ID:78384671)		Rumex dentatus	−9.9
Luteolin 7-rutinoside(PubChem ID:44258082)		Marrubium vulgare	−9.8
Kaempferol 7-(6″-galloylglucoside)(PubChem ID:74978085)		Acacia	−9.3

In order to gain more insight into the mechanism of interaction of these compounds with SARS-CoV-2 RdRp, a molecular docking was conducted using Autodock vina. Initially, Remdesivir which is the only approved drug to treat patients with severe COVID-19 was used to validate the docking approach (FDA, 12 April 2020). Remdesivir is a nucleotide analogue prodrug of (GS-441524) that converted into the active triphosphate form inside the cell. Recently, the structure of RdRp in complex with Remdesivir was solved (PDB: 7BV2). Remdesivir was found to interact covalently in its monophosphate form with the primer strand at the +1 position within the active site of SARS-CoV-2 RdRp (Yin et al., 2020). Fig. 2 a showed that the redocked Remdesivir monophosphate pose binds to the active site within RdRp as same as the co-crystalized ligand with RMSD = 2.969 Å. This result validates the ability of our docking protocol to predict the bioactive conformation of inhibitors. In addition, Remdesivir monophosphate exhibited lower binding free energy (−7.1 kcal/mol) in comparison to candidate compounds. This implies a strong possibility for these ligands to bind RdRp with higher affinity (Table 1). The molecular docking results of identified ligands illustrated that all compounds adapted binding modes similar to Remdesivir monophosphate within the active-site of RdRp (Fig. 2A). The chrysophanol 8-(6-galloylglucoside) was found to interact with the active site through six H-bonds (Tyr619, Asn691, Lys621, Asp623, Asp760 and Asp761) and one hydrophobic interaction (Pro621) (Fig. 2b). In the other hand, luteolin 7-rutinoside was involved in nine H-bonds with the active site (Arg553, Lys545, Arg555, Asp623, Asn691, Thr687, Asp760 and Ser759) (Fig. 2c). Fig. 2d, showed that kaempferol 7-(6″-galloylglucoside) formed seven H-bonds with RdRp (Arg553, Thr556, Tyr619, Trp617, Lys621, Asp760 and Asp761). Notably, the catalytic site of RdRp is composed of seven conserved motifs from A to G. Similar to Remdesivir, candidate compounds interact via H-bonds with several residues including the conserved Ser759-Asp60-Asp761 sequence within motif C of the RdRp polymerase domain (Fig. 2B) (Gao et al., 2020). These results suggested that the selected hits might act as non-nucleotide anti-polymerase agents.

Fig. 2

Binding mode and molecular interaction of compounds with RdRp. (A) Ribbon representation of RdRp (PDB: 7BV2) shown the binding mode of Remdesivir monophosphate (Liu et al., 2020), Chrysophanol8-(6-galloylglucoside) (yellow), Luteolin 7-rutinoside (orange) and Kaempferol 7-(6″-galloylglucoside) (blue). (B) The superimposition of co-crystalized (Liu et al., 2020) and docked (yellow) Remdesivir monophosphate (a). 2D ligand–protein interaction of (b) Chrysophanol8-(6-galloylglucoside), (c) Luteolin 7-rutinoside (Liu et al., 2020) and (d) Kaempferol 7-(6″-galloylglucoside).

In order to assess RdRp association and to determine the dynamic stability of screened compounds, all atom MD simulation at 20 ns were carried out by GROMACS. Analysis of RMSD from the starting structure revealed that most of complexes reach equilibrium within 20 ns (Fig. 3 A). This points to the high stability of protein–ligand complexes. Some fluctuation was observed within the first 10 ns in all systems. For the structure of RdRp bound to Remdesivir monophosphate and Chrysophanol8-(6-galloylglucoside), the RMSD values stabilized around ~0.32 nm. On the other hand, the RMSD values for the RdRp complex with Luteolin 7-rutinoside and Kaempferol 7-(6″-galloylglucoside), were equilibrated at around ~0.24 and ~0.20 nm, respectively. These results suggested that the binding of Remdesiver monophosphate and Chrysophanol8-(6-galloylglucoside) to RdRp may induce conformational changes (Fig. 3A). In consistent with this, the analysis of RMSF vs RdRp residue number showed that Remdesivir monophosphate and Chrysophanol8-(6-galloylglucoside) complexes showed higher oscillations in backbone residues in comparison to Luteolin 7-rutinoside and Kaempferol 7-(6″-galloylglucoside) systems (Fig. 3B).

Fig. 3

(A) RMSD of backbone atoms (C, Cα, and N) for RdRp-Ligand complex systems. (B) RMSF of backbone atoms for RdRp-Ligand complex systems. (C) Radius of gyration of backbone atoms (D) Number of H-bonds over 20 ns simulation time.

The RMSF values were high in some residues located in the flexible loop as well as alpha helices regions such as Thr261, Asp336, Glu431 and Asp824. Radius of gyrate (Rg) is another parameter to assess the compactness changes of a ligand–protein complex. Remdesivir monophosphate displayed a high Rg fluctuation in comparison to other compounds. This is consistent with the docking results of Remdesiver monophosphate that showed the lowest binding free energy. In contrast, no significant changes were observed in Rg values of Chrysophanol8-(6-galloylglucoside) and Kaempferol 7-(6″-galloylglucoside) indicating sustained compactness of both systems throughout the simulation (Fig. 3C). The Luteolin 7-rutinoside complex was the most stable and compact system with Rg value of ~2.94 nm. Next, the formation and stability of H-bonds were investigated over the simulation time (Fig. 3D). The H-bond properties are essential parameter in drug design due to their role in drug specificity, metabolism and absorption (Bitencourt-Ferreira et al., 2019). The results illustrated that Remdesiver monophosphate, Chrysophanol8-(6-galloylglucoside), Luteolin 7-rutinoside and Kaempferol 7-(6″-galloylglucoside) could form up to 5, 10, 11 and 10H-bonds, respectively. Therefore, the stability of complexes was maintained by H-bonds formation with active site residues. According to the docking and MD simulation analyses, Luteolin 7-rutinoside showed the best affinity towards RdRp in comparison to other compounds. It showed a high docking score (−9.8 kcal/mol) and was able to form stable H-bonds with the essential Asp760 and Asp761 of SDD sequence of polymerase domain. Moreover, the results from RMSD, RMSF and Rg indicated that the Luteolin 7-rutinoside complex was stable during the simulation period. Therefore, it can be concluded that Luteolin 7-rutinoside is the most preferable compound for the inhibition of RdRp for further in vitro and in vivo studies.

3-Chymotrypsin-like protease (3CLpro)

The 3CLpro, also known as the main protease (Mpro) is a cysteine protease. Its proteolytic activity is mediated by His41-Cys148 catalytic dyad (Anand et al., 2003, Needle et al., 2015). The structure based virtual screening of our phytochemical database led to identification of three compounds namely, 3,4,5-Tri-O-galloylquinic acid, chrysophanol 8-(6-galloylglucoside) and iso-mulbel-rochromene with docking scores of −9.7, −9.7- and −9.6 kcal/mol, respectively (Table 2 ). Interestingly, Chrysophanol 8-(6-galloylglucoside) was also identified as a potential RdRp inhibitor indicating that it might have a dual inhibitory mechanism.

Table 2

Potential inhibitors of 3CLpro from in-house Saudi medicinal database.

Name/ID	Chemical structure	Natural source	Docking score (kcal/mol)
X77(PubChem ID: 145998279)		–	−8.6
3,4,5-Tri-O-galloylquinic acid(PubChem ID: 127406)		Guiera senegalensis	−9.7
Chrysophanol 8-(6-galloylglucoside)(PubChem ID: 78384671)		Rumex dentatus	−9.7
Iso-mulbel-rochromene(PubChem ID: 196583)		Chenopodium ambrosioides, Morus alba L	−9.6

To predict these inhibitors conformation and orientation inside the active site of 3CLpro, a molecular docking was carried out. All selected compounds were found to fit into the active site with similar binding mode (Fig. 4 A). The co-crystalized non-covalent ligand, X77, was re-docked to assess the capability of our docking protocol to re-produce the bioactive conformation. As shown in Fig. 4a, the docked pose adapted similar binding mode as the co-crystalized ligand with RMSD = 0.789 Å indicating the robustness of our docking procedure. 3,4,5-Tri-O-galloylquinic acid binds with the active site of 3CLpro through several H-bonds (Met49, Leu141, Ser144, Glu166, His163 and Asp187) (Fig. 4b). The interaction between chrysophanol 8-(6-galloylglucoside) and 3CLpro was essentially mediated through H-bonds with four residues (Glu166, Arg188, Thr190 and Gln192) (Fig. 4c). Iso-mulbel-rochromene formed six H-bonds with the active site of 3CLpro (Thr25, Leu141, Ser144, His163, Arg188 and Gln189). These compounds block most of the essential residues involved in substrate recognition through H-bonds interaction indicating that these compounds may have potential to inhibit SARS-CoV-2 3CLpro activity.

Fig. 4

Binding mode and molecular interaction of compounds with 3CLpro. (A) Ribbon representation of 3CLpro (PDB: 6w63) shown the binding mode of co-crystalized (Pettersen et al., 2004) and docked (pink) X77, 3,4,5-Tri-O-galloylquinic acid (orange), chrysophanol8-(6-galloylglucoside) (cyan), iso-mulbel-rochromene (blue). (B) The superimposition of co-crystalized (Pettersen et al., 2004) and docked (pink) X77 (a). His41-Cys145 catalytic dyad is shown in cyan sticks. 2D ligand–protein interaction of (b) 3,4,5-Tri-O-galloylquinic acid, (c) chrysophanol8-(6-galloylglucoside) and (d) iso-mulbel-rochromene with 3CLpro.

Next, we characterized the stability of dynamic behavior and binding conformation of 3CLpro in complex with co-crystalized ligand, X77 as well as the identified compounds by 20 ns all-atom MD simulations. Fig. 5 A illustrated the plot of RMSD (nm) vs. time (ns) for apo-3CLpro and three 3CLpro-ligand complexes (3,4,5-Tri-O-galloylquinic acid, chrysophanol8-(6-galloylglucoside) and iso-mulbel-rochromene). All complexes attained the equilibration state during the 20 ns MD simulation period. The 3CLpro-X77 complex showed higher RMSD value as compared to identified ligand complexes during the last 10 ns, whereas 3,4,5-Tri-O-galloylquinic acid showed the least RMSD value, i.e., a greater dynamic stability compared to the other compounds. The RMSF values showed normal fluctuation with all ligand–protein complexes with respect to the RMSF value of co-crystalized ligand (Fig. 5B). The Rg profiles were quite similar for all ligand-complexes as compared to X77-3CLpro complex. Thus, protein-ligands complexes were stable complexes and well folded (Fig. 5C).

Fig. 5

(A) RMSD of backbone atoms (C, Cα, and N) for 3CLpro-Ligand complex systems. (B) RMSF of backbone atoms for 3CLpro-Ligand complex systems. (C) Radius of gyration (Rg) of backbone atoms (D) Number of H-bonds over 20 ns simulation time.

The number of H-bonds for ligand–protein complexes was calculated over 20 ns simulation time. Fig. 5D showed that 3,4,5-Tri-O-galloylquinic acid, chrysophanol8-(6-galloylglucoside) and iso-mulbel-rochromene had an average of eight, eleven and seven H-bonds, respectively, whereas an average of four H-bonds are formed in the X77-3CLpro complex. This explains the high binding energy scores of identified ligands compared to co-crystallized ligand, X77. These findings offer an extensive opportunity for these natural compounds, especially 3,4,5-Tri-O-galloylquinic acid, to contradict 3CLpro for the treatment of COVID-19.

Papain-like protease (PLpro)

The PLpro is another key cysteine protease enzyme for coronaviruses life cycle. In addition to its function in virus replication, it is involved in facilitating the virus entry by antagonizing the host cell immunity (Báez-Santos et al., 2015, Chen et al., 2020, Lei et al., 2018, Liu et al., 2020). The virtual screening against PLpro revealed that withanolide A, isocodonocarpine and calonysterone exhibited the highest binding affinity scores of −7.4, −7.0 and −6.9 kcal/mol, respectively (Table 3 ).

Table 3

Potential inhibitors of PLpro from in-house Saudi medicinal database.

Name/ID	Chemical structure	Natural source	Docking score (kcal/mol)
GRL0617(PubChem ID: 24941262)		–	−6.5
Withanolide A(PubChem ID: 11294368)		Datura innoxia	−7.4
Isocodonocarpine		Capparis decidua	−7.0
Calonysterone(PubChem ID: 101281312)		Senna obtusifolia	−6.9

Next, we performed a molecular docking to determine the molecular interactions of these inhibitors with PLpro. At first, the docking protocol was validated by docking a known SARS-CoV inhibitor; GRL0617 (PDB ID: 3E9S) into the same site within SRS-CoV-2 PLpro (PDB ID: 6W9C) (Ratia et al., 2008). This site is a conserved in PLpro for both SARS viruses. It is located in a proximity to the catalytic site which is composed of catalytic triad residues (Cys112-His273-Asp287) of PLpro (Fig. 6 A). Notably, the initial virtual screening was carried out against this pocket. Fig. 6A showed that GRL0617 as well as all identified ligands fit into the binding pocket with similar binding modes. A closer view into the binding of GRL0617 to PLpro of SRAS-CoV-2 showed that it adapted a similar binding mode and conformational arrangement for SARS-CoV (PDB ID: 3E9S), which demonstrated the validity of docking approach (Fig. 6a). Withanolide A interacted with the active site of PLpro through three H-bonds with Lys157, Asp164 and Glu167 (Fig. 6b). In the other hand, the cyclic compound, isocodonocarpine formed hydrophobic interactions with Tyr264 and Tyr268 and π-π stacking with Tyr264 (Fig. 6c). The calonysterone’s binding was found to involve in three H-bonds with Glu161, Asp164 and Tyr246 and one hydrophobic interaction with Pro248 (Fig. 6d).

Fig. 6

Binding mode and molecular interaction of compounds with PLpro. (A) Ribbon representation of PLpro (PDB: 6w9c) showing the binding mode of docked GRL0617 (Pettersen et al., 2004), withanolide A (Liu et al., 2020), isocodonocarpine (blue) and calonysterone (orange). The active site and catalytic triad Cys112-His273-Asp287 is shown in red sticks. (B) The superimposition of binding mode of docked GRL0617 (Pettersen et al., 2004) to SRAR-CoV-2 (cyan) (PDB: 6w9c) with the binding mode of co-crystalized GRL0617 (yellow) in complex with SARS-CoV (gold) (PDB: 3E9S). (a). 2D ligand–protein interaction of (b) withanolide A (c) isocodonocarpine and (d) calonysterone with PLpro.

Similarly, to ensure the stability of ligands within the active site of PLpro, 20 ns all-atom MD simulation using GROMAS 2018 was performed. The RMSD results of withanolide A-PLpro, isocodonocarpine-PLpro and calonysterone-PLpro complexes reached a stable state throughout the simulation period. The fluctuation in the RMSD values of withanolide A- and isocodonocarpine-PLpro complexes were less compared to the apo-protein and calonysterone-PLpro complex (Fig. 7 A). The RMSF values for isocodonocarpine showed less RMS fluctuation than other complexes, with respect to the RMSF values of apo-protein (Fig. 7B).

Fig. 7

(A) RMSD of backbone atoms (C, Cα, and N) for PLpro-Ligand complex systems. (B) RMSF of backbone atoms for PLpro-Ligand complex systems. (C) Radius of gyration (Rg) of backbone atoms (D) Number of H-bonds over 20 ns simulation time.

Fig. 7C illustrated that all ligand–protein complexes showed similar Rg behaviour, with respect to the apo-protein suggesting stable ligand–protein complexes in general. The H-bonds were also measured for each complex (Fig. 7D). Withanolide A, isocodonocarpine and calonysterone showed an average of two, three and five H-bonds. Although, isocodonocarpine did not show any H-bonds by molecular docking, three H-bonds were predicted by MD simulation. In addition, the number of H-bonds formed by withanolide A were almost stable, while the number of H-bonds formed by calonysterone was decreased over the MD simulation. This could explain the low docking score of calonysterone. Overall, Withanolide A showed the best results with respect to docking score as well as the MD simulation analyses suggesting that it might be a potential PLpro inhibitor candidate to combat COVID-19 infection.

Conclusion

In the presented study, we employed computational screening approach for the discovery of potential anti-COVID-19 against the novel coronavirus SARS-CoV-2 infection. A total of ~1000 phytochemicals derived from 60 Saudi medicinal plants with recognised antiviral activities have been investigated and lead to the identification of nine potential inhibitors for three major druggable targets of SARS-CoV-2. The obtained candidates from the virtual screening were subjected to MD simulation to confirm their binding and dynamic stability at the target site. It has been demonstrated that chrysophanol 8-(6-galloylglucoside), luteolin 7-rutinoside and kaempferol 7-(6″-galloylglucoside) have the potential for further optimisation as an RdRp inhibitors. Moreover, 3,4,5-tri-O-galloylquinic acid, chrysophanol 8-(6-galloylglucoside) and iso-mulbel-rochromene showed potential to efficiently bind and inhibit the viral 3CLpro. Withanolide A, isocodonocarpine and calonysterone have also shown promising results in targeting the PLpro. Future in vitro and in vivo evaluations of these compounds can confirm their antiviral activity and subsequently their potential to serve as drug candidates for further designing of target-specific anti-COVID-19 agents.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.