Potential role of medicinal plants and their constituents in the mitigation of SARS-CoV-2: identifying related therapeutic targets using network pharmacology and molecular docking analyses.

Since the outbreak of Coronavirus disease (COVID-19) caused by SARS-CoV-2 in December 2019, there has been no vaccine or specific antiviral medication for treatment of the infection where supportive care and prevention of complications is the current management strategy. In this work, the potential use of medicinal plants and more than 16 500 of their constituents was investigated within two suggested therapeutic strategies in the fight against SARS-CoV-2 including prevention of SARS-CoV-2 RNA synthesis and replication, through targeting vital proteins and enzymes as well as modulation of the host's immunity through production of virulence factors. Molecular docking studies on the viral enzymes 3Clpro, PLpro and RdRp suggested rocymosin B, verbascoside, rutin, caftaric acid, luteolin 7-rutinoside, fenugreekine and cyanidin 3-(6''-malonylglucoside) as promising molecules for further drug development. Meanwhile, the medicinal plants Glycyrrhiza glabra, Hibiscus sabdariffa, Cichorium intybus, Chrysanthemum coronarium, Nigella sativa, Anastatica hierochuntica, Euphorbia species, Psidium guajava and Epilobium hirsutum were enriched in compounds with the multi-targets PTGS2, IL2, IL1b, VCAM1 and TNF such as quercetin, ursolic acid, kaempferol, isorhamnetin, luteolin, glycerrhizin and apigenin. Enriched pathways of the molecular targets included cytokine-cytokine receptor interaction, TNF signaling pathway, NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, NF-kappa B signaling pathway and JAK-STAT3 signaling pathway which are all closely related to inflammatory, innate and adaptive immune responses. The present study identified natural compounds targeting SARS-CoV-2 for further in vitro and in vivo studies and emphasizes the potential role of medicinal plants in the mitigation of SARS-CoV-2. This journal is © The Royal Society of Chemistry.

Introduction

Coronaviruses (CoV) are a family of enveloped viruses containing a large positive-sense single-stranded RNA. They are distributed in a wide range of animals and can infect humans, primarily afflicting the pulmonary system and gastrointestinal tract. Four genera of CoV exist: alpha-, beta, gamma- and delta coronavirus. Alpha- and beta-coronavirus stem from mammals, particularly bats, while gamma- and delta-coronaviruses are of avian origin.[1]

In recent years, the world experienced two pandemics due to CoV. These outbreaks include the severe acute respiratory syndrome coronavirus (SARS-CoV) and the Middle East respiratory syndrome coronavirus (MERS-CoV). During late December of 2019, a cluster of hospital admitted patients were diagnosed with severe pneumonia of unknown etiology. This cluster was epidemiologically linked to a local seafood market in Wuhan, Hubei Province, People's Republic of China.[2,3] These patients were found to be infected with the novel coronavirus strain, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2 – previously named 2019-nCoV) which belongs to the beta-coronavirus genera.[4]

Infection with SARS-CoV-2 can be asymptomatic or gives rise to a symptomatic disease termed Coronavirus Disease 2019 (COVID-19) which varies in severity and may lead to fatality.[5] The median incubation period is reported to be 5.1 days and almost all infected symptomatic patients show symptoms in no more than 12 days of infection.[6] In cases of mortality, the number of days from symptom onset to death was reported to range from 6 to 41 days with an average of 14 days.[7] The most commonly occurring symptoms of COVID-19 include fever, dry cough, and fatigue.[8,9] Initially, some patients may present with gastrointestinal upset in the form of diarrhea and nausea prior to onset of fever.[9] Other symptoms include dyspnoea, sputum production, headache, haemptyosis[10] and lymphopenia.[8] The most common complication during hospitalization is pneumonia, preceding acute respiratory distress syndrome (ARDS).[8]

Another side of SARS-CoV-2 infection is the cytokine storm syndrome which is brought upon the body as a result of the immune system's excessive response. The culprits, cytokines and chemokines are expressed and secreted upon pathogen encounter and although their importance in the body's immunity is well-documented, dysregulation and abundance in their production does more harm than good.[11] COVID-19 patients were found to have an elevated level of interleukin 1B (IL-1B), IL-6, IL-10, interferon gamma (IFN-γ), tumor necrosis factor α (TNF-α), interferon gamma-induced protein (IP10), granulocyte colony-stimulating factor (GCSF) and MCP1 (monocyte chemoattractant protein-1). In addition, ICU admitted patients showed higher levels of plasma cytokines compared to their non-ICU admitted counterparts; indicating that the severity of COVID-19 illness is amplified by cytokine levels.[10,12] Moreover, serum levels of IL-6, IL-10 and TNF-α in non-ICU admitted patients showed negative correlation with T-cell count. This infers that higher levels of these cytokines elicit a negative regulatory response in the production of T-cell or their survival.[12] Hence, inhibition of specific pro-inflammatory cytokines may lessen the hyper-inflammation involved in COVID-19 pathogenesis.[13]

In patients suffering from COVID-19, SARS-CoV-2 binds to epithelial cells present in the alveoli, thus invoking and activating the body's innate and adaptive immune system. This leads to the production and release of cytokines in order to further stimulate the immune system. As a result, an increase in blood vessel permeability occurs and subsequently, blood cells and fluids seep into the alveoli and this results in dyspnea.[14,15] If uncontrolled, a ramification of continued hyper-inflammation and diffuse damage of the alveoli is ARDS.[16]

Currently, there are no vaccines or specific antiviral medication to treat SARS-CoV-2 infections.[10] Therefore, supportive care and prevention of complications is an important management strategy to minimize harm. Suggested therapeutic strategies in the fight against SARS-CoV-2 can be split into 4 different categories: (1) prevention of SARS-CoV-2 RNA synthesis and therefore replication, through action on vital proteins and enzymes, (2) impeding virus-cell receptor binding or suppression of virus self-assembly, (3) stimulation of the host's immunity through production of virulence factors and (4) blocking the entry of the virus into host cells through action on the host's enzymes or cellular entry receptor.[17]

Nonstructural proteins (Nsp) are functional proteins of high importance to CoV. They take part in the process of virus replication and host infection through their role in RNA transcription and translation and the subsequent processing and modification of synthesized proteins.[17] Nsp are produced through proteolytic cleavage of replicase polyprotein 1a (pp1a) and replicase polyprotein 1ab (pp1ab) by the action of viral papain-like protease (PLpro) on the N-termini resulting in 3 products[18] and 3-chymotrypsin-likeprotease (3CLpro, also named main protease) at central region and c-termini which gives rise to a minimum of 12 products, 3CLpro included.[19] From the aforementioned cleavage, a total of 16 Nsp are produced which are named Nsp1–Nsp16.[20] Nsp12 is an RNA-dependent RNA polymerase (RdRp) enzyme that is vital in the lifecycle of CoV as it is responsible for the replication–transcription complex.[21,22] However, despite the significant role that RdRp plays in RNA viruses, there is little data on RdRp inhibitors in any CoV. Therefore, it is clear that these key viral enzymes are of great importance to CoV. And owing to their essential active sites, it is hypothesized that small inhibitory molecules are of high therapeutic potential. Other targets include viral spike surface glycoprotein (S) which has a critical role in virus-cell receptor binding. However, this is a target for larger molecules such as lipopeptide EK1C4 as reported by Xia, S. et al.[23]

For centuries, natural products have been a rich source of medicinally active constituents that have had a pivotal role in treating and preventing an infinite number of diseases. Natural products are a cornerstone of Traditional Chinese Medicine (TCM), which encompasses a tremendous wealth of ancient knowledge. The importance of TCM in emerging pandemics was made evident when it was effectively partnered with conventional medicine in the management of the SARS outbreak in 2003.[24] Therefore, due to their notable contributions in the SARS-CoV pandemic, the National Health Commission of the People's Republic of China employed TCM in the COVID-19 treatment guidelines.[25] At a press conference held on February 17, an official reported that 60 107 COVID-19 cases (85.2% of total cases at the time) were treated with TCM.[26]

Additionally, in an investigation covering 23 Chinese provinces, it was reported that the most frequently prescribed Chinese herbs used in the prevention of COVID-19 were (in descending order): Radix astragali, Glycyrrhizae Radix Et Rhizome, Radix Saposhnikoviae, Rhizoma Atractylodis macrocephalae, Lonicerae japonicae Flos, Fructus forsythia, Atractylodis Rhizoma, Radix Platycodonis, Pogostemonis Herba, Cyrtomium fortunei J. Sm.[27]

Within this context, the potential of the use of medicinal plants in the mitigation of the novel SARS-CoV-2 infection has been investigated. A database comprising more than 16 000 compounds was compiled and docked against the crucial viral proteins; 3-chymotypsin protease (3CLpro), papain-like protease (PLpro) and RNA-dependent RNA polymerase (RdRp) to find potential inhibitors for these enzymes. Furthermore, network pharmacology analysis of all the compounds in the database was attempted to speculate those that can target the inflammatory and immunity-related pathways at the molecular level and revealing those plants with potential multi-components and multi-targets that may help regulate the body functions and possibly will play a therapeutic role in the mitigation of the disease.

Materials and methods

Molecular docking studies

The Protein Data Bank (PDB) was utilized to retrieve the crystal structure of the three SARS-CoV-2 viral proteins; main proteinase or chymotrypsin-like protease (3CLpro, PDB ID: 5R7Y), papain-like protease (PLpro, PDB ID: 6W9C) and RNA-Dependent RNA Polymerase (RdRp, PDB ID: 6M71). The proteins structures were prepared by means of protein preparation module of Schrodinger's Maestro Molecular modeling suit (Schrödinger Release 2018-1). The proteins were first preprocessed by assigning bond orders and hydrogens in addition to creation of zero order bonds to metals and disulphide bonds, all the water molecules were deleted beyond 5 Å from the active site, hydrogen bonds were assigned using PROPKA at PH = 7 followed by energy minimization using OPLS 3 force field until the relative mean standard deviation (RMSD) of the minimized structure compared to the crystal structure exceeded 0.30 Å.[28] Location of the binding site for the docking experiments was determined using the receptor grid generation module and the boxes enclosing the centroids of co-crystallized ligands were set as the grids.

The 3D structures of the compounds were imported as SDF file into the LigPrep module of the Maestro molecular modeling package to deliver low energy structures of compounds. Ionization states were set to generate all possible states at pH 7. Molecular docking simulations were performed using Glide docking program of Maestro molecular modeling package implementing SP- and XP-Glide module. Compound–target interactions including hydrogen bond, ion pair interactions, hydrophobic interactions and the binding modes were analyzed in Maestro interface. The docking scores of the top 50 interacting compounds were utilized for construction of another compound–target network using Cytoscape 3.5.1 (http://www.cytoscape.org/) software.

Study design of network pharmacology analysis

The 2D structures of 16 500 compounds from medicinal plants were imported from different resources including NANPDB (http://african-compounds.org/nanpdb) public database, Dictionary of Natural Products (DNP) and two in-house database.[28,29] The chemical structures of the compounds were confirmed through ChEMBL (https://www.ebi.ac.uk/chembl/) and PubChem (https://pubchem.ncbi.nlm.nih.gov/). Schrodinger software (2018-1) was used to convert the 2D chemical structures of the compounds to SMILES format. Qikprop software (Schrodinger suite 2018-1) was utilized for drug-likeness study to filter out the compounds by calculation of absorption, distribution, metabolism, and excretion (ADME) criteria and Lipinski's rule of five.

Identifying potential targets/genes of the filtered chemical compounds was achieved through the public database STITCH DB (http://stitch.embl.de/, ver. 5.0) and Swiss Target Predication (http://www.swisstargetprediction.ch) with the ‘Homo sapiens’, while STRING database (https://string-db.org/) was used to draw a protein–protein interaction network (PPI network). The identified targets were further filtered by retaining only those correlated to inflammation or immunity through searches in the UNIPROT and Therapeutic Target Database (TTD).

Information about gene ontology including the highly-associated pathways, biological processes and functional annotation of the target protein was retrieved from Database for Annotation, Visualization and Integrated Discovery (DAVID) ver. 6.8 (https://david.ncifcrf.gov/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (http://www.genome.jp/kegg/pathway.html). Pathways with P-values < 0.05 were only retained.

The obtained results were employed for construction of Compound–Target, Compound–Target–Pathway and Compound–Target–Pathway-Plant visual networks using Cytoscape software. The network parameters were calculated using the network analyzer plug-in in Cytoscape. The importance of nodes in each network was judged through their combined score of interactions.

Results and discussion

Identifying potential medicinal plants constituents against viral targets

A medicinal plant/natural products library comprised of 10 823 phytoconstituents was screened in silico via molecular docking against three proven protein targets in SARS-CoV-2 striving for potent anti-COVID-19 natural compounds. 3-Chymotrypsin-like protease (3CLpro), papain-like protease (PLpro) and RNA-dependent RNA polymerase (RdRp) are elemental for viral replication and therefore, impending targets for anti-coronaviruses screening. Lopinavir, nelfinavir and remdesivir are commonly used antiviral drugs acting on viral proteins that are presently available in the market and undergoing clinical trials for repurposing in treatment of SARS-CoV-2 infection. These drugs were also analyzed against the three selected key protein targets. Meanwhile, other drugs like hydroxychloroquine and dexamathasone don't act on the viral proteins but rather on the immune system reaction caused by the viral infection so they were not involved in the docking study.

Molecular docking against 3-chymotrypsin-like protease (3CLpro)

The main protease (Mp) or 3-chymotrypsin-like protease (3CLpro) or Nsp5 has a critical role in cleavage of the polyprotein at eleven distinct sites to generate various Nsps that are essential for viral replication.[30] The maturation of Nsps, which is crucial in the virus life cycle is directly mediated via 3CLpro. The detailed inspection on catalytic mechanism of 3CLpro makes it an attractive target for anti-COVID-19 drug development.[31]

Virtual screening of our database using Glide Xp protocol against 3CLpro crystalline structure (5R7Y) deposited in March in complex with co-crystallized ligand disclosed that rutin from Anastatica hierochuntica, Glycyrrhiza glabra and Allium myrianthum followed by rocymosin b from Glycyrrhiza glabra and subsequently verbascoside from Cichorium intybus, Olea europaea and Marrubium vulgare were the top three scoring compounds (Table 1) which endorsed us to investigate their binding modes. The compounds bound to the Cys–His catalytic dyad (Cys-145 and His-41)[32] along with the other residues, with notable docking scores (Table 1). It was discerned that the active site of the main protease in which the three top scoring compounds bind enclosed hydrophobic residues such as Met49, Leu141, Cys145, Met165, Leu167 and Pro168 in addition to amino acids polar contributions such as Thr25, Thr26, His41, Asn142, Ser144, His163 and Hie164. Hydrogen bonding with Leu141 and Glu166 also stabilized their conformations (Fig. 1).

Docking scores of the top 20 hits identified in addition to lopinavir, nelfinavir and remedesivir against the three viral proteins studied

	3CLpro			PLpro			RdRp
	Compound	Docking score	Plant source	Compound	Docking score	Plant source	Compound	Docking score	Plant source
1		−12.632	Anastatica hierochuntica, Allium myrianthum, Glycyrrhiza glabra		−14.041	Cichorium intybus, Marrubium vulgare, Olea europaea		−11.541	Cichorium intybus
2		−11.844	Glycyrrhiza glabra		−13.883	Cynara scolymus		−10.664	Marrubium vulgare, Olea europaea, Cichorium intybus
3		−11.721	Cichorium intybus, Marrubium vulgare, Olea europaea		−11.148	Marrubium vulgare, Olea europaea, Cichorium intybus		−9.894	Trigonella foenum-graecum
4		−11.624	Cichorium intybus, Hibiscus sabdariffa		−10.832	Glycyrrhiza glabra		−9.797	Tribulus terrestris
5		−11.412	Glycyrrhiza glabra		−10.674	Ephedra alata		−9.785	Eruca sativa
6		−11.409	Sonchus macrocarpus		−10.658	Salvia triloba		−9.754	Hibiscus sabdariffa
7		−11.363	Glycyrrhiza glabra		−10.593	Glycyrrhiza glabra		−9.747	Matricaria chamomilla
8		−11.119	Terminalia species, Tamarix nilotica, Quercus suber		−10.504	Centaurea furfuracea		−9.732	Glycyrrhiza glabra
9		−11.061	Hibiscus sabdariffa		−10.49	Artemisia judaica		−9.611	Glycyrrhiza glabra
10		−11.015	Artemisia herba-alba		−10.402	Matricaria chamomilla		−9.544	Medicago sativa
11		−11.015	Epilobium hirsutum		−10.338	Phyllanthus emblica, Globularia alypum		−9.527	Tribulus terrestris
12		−10.892	Schinus molle		−10.255	Phyllanthus emblica		−9.509	Cleome species
13		−10.88	Lepidium sativum, Cleome species, Tribulus terrestris		−10.201	Glycyrrhiza glabra		−9.493	Trigonella foenum-graecum
14		−10.831	Schinus molle		−10.182	Tetraclinis articulata		−9.485	Fraxinus oxycarba
15		−10.811	Crataegus sinaica		−10.18	Olea europaea		−9.484	Glycyrrhiza glabra
16		−10.762	Olea europaea		−10.162	Cynara scolymus		−9.403	Artemisia judaica
17		−10.695	Olea europaea		−10.041	Podocarpus gracilior		−9.361	Glycyrrhiza glabra
18		−10.684	Cicer arietinum		−10.035	Glycyrrhiza glabra		−9.31	Scrophularia saharae
19		−10.658	Hibiscus sabdariffa		−10.027	Anastatica hierochuntica		−9.298	Salvia triloba
20		−10.624	Artemisia herba-alba, Artemisia judaica		−10.027	Glycyrrhiza glabra		−9.278	Tribulus terrestris
		−10.195		Lopinavir	−8.625		Lopinavir	−5.94
		−8.297		Nelfinavir	−4.324		Nelfinavir	−7.594
		−7.346		Remdesivir	−3.947		Remdesivir	−8.526

Fig. 1

3D (to the left) and 2D (to the right) interaction diagrams of (A) rutin, (B) rocymosin b and (C) verbascoside with the crystal structure of 3CLpro (PDB ID: 5R7Y).

Molecular docking against papain-like protease (PLpro)

PLpro is the proteolytic enzyme that cleaves the N-terminus of viral replicase poly-protein at three sites to release mature Nsp1, 2 and 3. This process is fundamental for correcting virus replication.[18] Moreover, PLpro also has a critical role in antagonizing the host's innate immunity.[33] As a result of its imperative function in the process of coronavirus replication and host infection, PLpro has emerged as a prevalent target for coronavirus inhibitors. However, no PLpro inhibitor has reached the FDA marketing approval till this moment.

After conducting the extra precision Glide docking protocol on our database and the three chosen antiviral repurposed drugs, the interactions between the three top scoring constituents and PLpro crystalline structure (PDB ID: 6W9C) released on first of April 2020, were analyzed (Fig. 2 and S2†). It was observed that the specific pocket in which the top scoring compounds bind contained non-polar amino acids like Leu162, Val202, Met206, Met208 and Ile222 in addition to amino acids polar contributions such as Ser170 and Asn156.

Fig. 2

3D (to the left) and 2D (to the right) interaction diagrams of (A) verbascoside, (B) luteolin-7-rutinoside and (C) caftaric acid with the crystal structure of PLpro (PDB ID: 6W9C).

Verbascoside from Cichorium intybus, Olea europaea and Marrubium vulgare which achieved the highest XP docking score (Table 1) was furthermore bound to PLpro active site through hydrogen bonding with Glu167, Tyr171, Met208, Tyr268, Tyr273 and Thr301 and additionally through charged negative interactions with Glh203 and Asp164 (Fig. 2A).

It was also discerned that verbascoside had bound to the two selected proteases PLpro and 3CLpro with high XP docking scores. These dual targets inhibition advocates verbascoside as a potential druggable anti-COVID-19.

The second active top scoring compound against PLpro; luteolin-7-rutinoside from Cynara scolymus (Table 1) interacted with numerous hydrogen bonds with the amino acids Asn156, Gly163, Ser170, Tyr268, Tyr273 and Thr301 as well as π–π stacking interaction with Arg166 and negatively charged interactions with Glh161, Asp164, Glh203 and Asp302 (Fig. 2B).

Analysis of binding modes of caftaric acid; the third top scoring compound, with PLpro revealed that besides the polar and non-polar amino acids contributions in the active site, caftaric acid exhibited four H-bonding with Leu162, Tyr207, Met208 and Lys232 (Fig. 2C).

Molecular docking against RNA-dependent RNA polymerase (RdRp)

RNA-dependent RNA polymerase (RdRp) or nonstructural protein 12 (Nsp12) is the crucial enzyme for coronavirus replication/transcription complex. It catalyzes the replication of RNA from an RNA template upon binding to cofactors. The RdRp domain of polymerase is located at the C-terminus and possesses a conserved Ser–Asp–Asp motif.[33] The crystal structure of SARS-Cov-2 RdRp in complex with cofactors (PDB ID: 6M71) has been recently released in RCSB PDB on the first of April 2020 which provided an excellent ground for structure-based drug discovery.[34] Targeted inhibition of RdRp could not cause host cells significant side effects or toxicity.[35] Fast drug repurposing attempts have recognized the antiviral remdesivir, a RdRp inhibitor, as a promising anti-COVID-19.[36] Clinical trials are currently continuing to investigate the full efficiency spectrum of remdesivir in COVID-19 patients.[37]

Molecular docking studies of the compiled database against RdRp crystal structure (PDB ID: 6M71) revealed that the glycoside “cyanidin 3-(6''-malonylglucoside)” present in Cichorium intybus showed the lowest binding energy with XPG score of −11.541 kcal mol−1. As shown in Fig. 3A, cyanidin 3-(6′′-malonylglucoside) fitted well into RdRp active site. Eight hydrogen bonds were predicted between Lys545, Arg553, Arg555, Tyr619, Asp623, Asp760 and the carbonyl and hydroxyl groups of the compound. In addition, charged negative interactions were also observed with Asp 618, 623, 760 and 761.

Fig. 3

3D (to the left) and 2D (to the right) interaction diagrams of (A) cyanidin 3-(6′′-malonylglucoside), (B) caftaric acid and (C) fenugreekine with the crystal structure of RdRp (PDB ID: 6M71).

The second top scoring compound “caftaric acid”, a constituent of Cichorium intybus, Olea europaea and Marrubium vulgare, showed a dual SARS-CoV-2 target inhibition; PLpro as previously described with RdRp. It possessed a docking score of −10.664 kcal mol−1. The binding modes of its interaction with RdRp are depicted in Fig. 3B, in which many of the hydrophobic amino acid residues such as Tyr455, Ala554, Val557 and Cys622 compose a relatively hydrophobic environment to accommodate the compound. Moreover, π–π stacking interactions with Arg553, Arg624, hydrogen bonding with Lys545, Arg553, Thr556, Asp623 and polar interaction with Ser682 stabilize its conformation within the active pocket.

Fenugreekine (from Trigonella foenum-graecum) was the third scoring compound with binding energy of −9.894 kcal mol−1. Analysis of binding modes (Fig. 3C) revealed numerous hydrogen bonds between fenugreekine and RdRp active site residues such as Trp617, Asp623, Asp760, Asp761, Lys551, Arg553, Arg555 and Thr556 in addition to Van der Waals forces with Val557, Cys622 and Trp619 residues. Polar interaction with Ser549, Thr556 and Ser682 as well as charged negative interactions with Asp542, Asp618, Asp623 and Asp761 were also perceived.

It is noteworthy that these natural top-scored phytoconstituents surpassed in their binding free energies (Table 1) the three selected antiviral drugs adopted from FDA drug repurposing criteria; remdesivir; the proven RdRp inhibitor which is currently in stage III clinical trial, and the anti-HIV drugs; lopinavir and nelfinavir. The 3D and 2D interaction diagrams of the three repurposed anti-COVID-19 drug candidates with RdRp, PLpro and 3CLpro crystalline structures are displayed in the ESI Fig. S1–S3.† The difference in stereochemistry observed in the docked poses is due to the ligand preparation module of Schrodinger suite which generates the chemical structure of the compound at physiological pH which may cause change in ionization or stereochemistry of the compound.

Our analyses revealed that theses top scoring compounds of natural compounds might serve as potential druggable anti- SARS-CoV-2 agents for further in vitro and in vivo studies and for drug repositioning, optimization and development to combat COVID-19.

Multi-target constituents and medicinal plants identification through docking network

The docking sores of the top 50 compounds against each individual antiviral target protein were exploited for construction of a visual network (Fig. 4). The nodes were sized according to their relative docking scores. As indicated from the size of its node shown in figure, 3CLpro enzyme was the protein target with the lowest docking scores and favorable binding energies followed by PLpro and finally RdRp. It can be observed that rocymosin B (a dihydrochalcone glycoside) was the only compound commonly observed in the top 20 scoring compounds of all three mutual antiviral targets; 3CLpro, PLpro and RdRp with docking scores of −11.844, −10.593 and −9.361, respectively. Meanwhile, verbascoside (acetoside, a caffeoyl phenylethanoid glycoside) exhibited remarkably low binding energy towards both 3-chymotrypsin-like and papain protease enzymes (Fig. 4) with docking scores of −11.721 and −14.041, respectively.

Fig. 4

Compound–viral target network colored according to number of common neighbors. Multi-target nodes are colored orange. Node size is proportional to the relative docking score of the compound in the network.

Investigating the plant sources of the top 50 scoring compounds in the docking study against the three antiviral protein targets revealed that the plants Glycyrrhiza glabra (liquorice), Cichorium intybus (chicory), Olea europaea (olive leaves), Marrubium vulgare (white horehound), Lepidium sativum (garden cress) and Artemisia judaica (Judean wormwood or Shih Balady) enclosed compounds that are mutually active against the three anti-viral targets (Fig. 5). The node size of Glycyrrhiza glabra (liquorice) indicated that it was enriched in the highest number of compounds possessing favorable binding poses in all three viral proteins, followed by Cichorium intybus, Olea europaea and Marrubium vulgare. Meanwhile Epilobium hirsutum (hairy willow-herb), Cicer arietinum (chickpea) and Centaurea incana (Kantarioun) possessed compounds with mutual targets against 3CLpro and PLpro. On the other hand, Hibiscus sabdariffa (Roselle), Tribulus terrestris, Trigonella foenum-graecum, Cleome species and Daucus carota possessed compounds that are mutually active against the two viral proteins 3CLpro and RdRp. Moreover, a dual PLpro, RdRp target binding was achieved by compounds from Matricaria chamomilla (German chamomile), Apium graveolens (celery), Medicago sativa (alfalfa) and Salvia triloba.

Fig. 5

Plant-viral target network colored according to number of common neighbors. Nodes are colored according to their topological coefficients; nodes with three targets are colored green and nodes with two targets are colored blue. Node size is proportional to the relative total docking scores of the compounds in the plants included the network.

Identifying potential targets of medicinal plants in inflammation and immune regulation pathways

Huge inflammatory response is observed in COVID-19 and in several clinical trials, it has been shown that combination of anti-inflammatory and antiviral drugs may be more effective than each one individually.[38] One of the main mechanisms involved in the SARS-CoV-2 induced-ARDS is the cytokine storm, the deadly uncontrolled systemic inflammatory response resulting from the release of large amounts of pro-inflammatory cytokines (IFN-α, IFN-γ, IL-1β, IL-6, IL-12, IL-18, IL-33, TNF-α, TGFβ, etc.). Therefore, arbitration of the immune evasion of SARS-CoV-2 is imperative in its treatment and specific drug development.[39-41] Targeting this approach by medicinal plants and their constituents and identifying their potential inflammation and immunity-correlated target proteins and pathways was the focus of this part of the study.

About 16 500 compounds from more than 1300 plant sources were compiled from different database including the North African plants database, Dictionary of Natural products and some in-house available databases. Duplicate compounds were removed and the resulting compounds were filtered using ADME screening based on their oral bioavailability (OB) and Lipinski's rule “rule of five”. Compounds having OB of less than 30% and disobeying two or more of the Lipinski's rule “rule of five” criteria were excluded from further analysis resulting in a database of 10 823 compounds.

Screening the compounds for possible target was performed through the STITCH 5.0 public databases and Swiss Target Prediction database, while the correlation of the identified targets to inflammation and immunity were retrieved from UniProt and TTD.

25 common proteins related to inflammatory and immunity pathways were identified as potential targets for 121 plants constituents (Table 2) through 1073 interactions in the constructed compound–target network with an average of 7.2 targets for each constituent (Fig. 6A). Protein correlation results indicated that the identified potential target proteins in the anti-inflammatory effects were correlated to and regulated each other (Fig. 6B).

Identified potential protein targets of the different compounds in medicinal plants and the involved pathways

Protein	Interacting compound (s)	Involved pathways
AKT1	Isorhamnetin, salidroside, apigenin, myricetin, genistein	VEGF signaling pathway, Toll-like receptor signaling pathway, B cell receptor signaling pathway, TNF signaling pathway, tuberculosis, influenza A
NOS2	Isorhamnetin	Pertussis, tuberculosis
MAPK8	Luteolin, quercetagetin, isorhamnetin	Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, TNF signaling pathway, pertussis, tuberculosis
PTGS2	Catechin, eupafolin, withaferin A, apigenin, medicocarpin, dl-liquiritigenin, licofuranocoumarin, astragalin, calycosin, formononetin, 6-prenylated eriodictyol	NF-kappa B signaling pathway, VEGF signaling pathway, TNF signaling pathway, prostaglandin-endoperoxide synthase activity
ALOX5	Protocatechuic acid, kaempferol, quercetin	Arachidonic acid metabolism
CASP3	Alternariol, apigenin, kaempferol, salidroside, gallic acid, isoquercitrin, myricetin, chrysin, hispolon	TNF signaling pathway, pertussis, tuberculosis
CASP8	Aloe-emodin	NOD-like receptor signaling pathway
CASP9	Aloe-emodin, alternariol	VEGF signaling pathway, tuberculosis, influenza A
FOS	Luteolin	Toll-like receptor signaling pathway, T cell receptor signaling pathway, B cell receptor signaling pathway, TNF signaling pathway, pertussis
HSP90AA1	Aloe-emodin	NOD-like receptor signaling pathway
IL1B	Quercetin, ursolic acid, isorhmnetin,	Cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, TNF signaling pathway, MAPK signaling pathway, IL-17 signaling pathway, Th17 cell differentiation, pertussis, influenza A, tuberculosis
IL2	Quercetin, ursolic acid	T cell receptor signaling pathway, Th17 cell differentiation, cytokine–cytokine receptor interaction
IL10	Salidroside	FoxO signaling pathway
STAT1	Acacetin, quercetin	Cytokine–cytokine receptor interaction
JUN	Acacetin, gallic acid, kaempferol, luteolin	Toll-like receptor signaling pathway, T cell receptor signaling pathway, B cell receptor signaling pathway, TNF signaling pathway
KDR	1-Methoxyphaseollidin, astragalin, gancaonin B, glycyrin, licochalcone G	Cytokine–cytokine receptor interaction, VEGF signaling pathway, MAPK signaling pathway, ras signaling pathway
PTGS2	Astragalin, calycosin, glyzaglabrin, isoliquiritigenin, isorhamnetin, kaempferol, naringenin, quercetin, ursolic acid	Arachidonic acid metabolism
TNF	Glycyrrhizin, kaempferol, naringin, quercetin, ursolic acid	Cytokine–cytokine receptor interaction, NF-kappa B signaling pathway, TNF signaling pathway, MAPK signaling pathway, IL-17 signaling pathway, HTLV-I infection, pertussis, human papillomavirus infection, T cell receptor signaling pathway, influenza A, tuberculosis
SMAD2	Luteolin	FoxO signaling pathway
VCAM1	Isoliquiritigenin, kaempferol, quercetin	NF-kappa B signaling pathway, TNF signaling pathway
VEGFA	Acacetin, hispidulin, quercetin, ursolic acid	NF-kappa B signaling pathway, VEGF signaling pathway, MAPK signaling pathway, human papillomavirus infection, PI3K-Akt signaling pathway
PIK3CG	Myricetin	TNF signaling pathway
NOS3	Genistein kaempferol, luteolin, luteolin-7-O-g, silibinin	VEGF signaling pathway, HIF-1 signaling pathway
GADPH	Koningic acid	HIF-1 signaling pathway
MMP9	Luteolin, gallic acid, hispolon	TNF signaling pathway

Fig. 6

(A) Compounds-target–pathway network. The orange circles represent the identified compounds, the green circles represent potential targets in inflammation and immunity. (B) Protein–protein interaction (PPI) network of identified targets colored according to degree of interaction.

The targets which possessed the highest combined score were presumed to be the most enriched in the network and were revealed to be PTGS2, IL2, IL1b, VCAM1 and TNF. Meanwhile, the constituents with the highest combined score were quercetin, ursolic acid, kaempferol, isorhamnetin, luteolin, glycerrhizin and apigenin (Fig. 7).

Fig. 7

(A) Analysis of the top 10 scoring plants merged plant-compound–target–pathway network. (B) The distributions % of the C–T interactions on the compounds identified. (C) The distributions % of the C–T interactions on the top 10 scoring medicinal plants in the database.

PTGS2 (prostaglandin-endoperoxide synthase 2) is a key enzyme in prostaglandin biosynthesis and is responsible for production of inflammatory prostaglandins. It has been recently hypothesized that the lipid mediator prostaglandin E2 may play a crucial role in SARS-CoV-2 infection progress and can be considered a pliable treatment intervention.[42]

IL2 (Interleukin-2 receptor subunit beta) and ILb (Interleukin-1 beta) are potent pro-inflammatory cytokines which are considered as the major endogenous pyrogen. They induce the synthesis of prostaglandins, neutrophil influx, T-cell and B-cell activation and cytokine production. Meanwhile, TNF (tumor necrosis factor) cytokine is mainly secreted by macrophages and is a potent pyrogen causing fever by direct action or by stimulation of interleukin-1 secretion and impairs regulatory T-cells function. VCAM1 (vascular cell adhesion protein 1) has an important role in leukocyte-endothelial cell adhesion and may play a pathophysiologic role both in immune responses and in leukocyte migration to sites of inflammation. The identified molecular targets are strongly related to the pathophysiology of SARS-CoV-2 infection where it has been shown that most patients with severe SARS-CoV-2 exhibit-noticeably high serum levels of the pro-inflammatory cytokines IL-6 and IL-1β, as well as IL-2, IL-8, IL-17, and TNF, which leads to what is known as the ‘cytokine storm’. It is postulated that small-molecule inhibitors of their downstream signaling components may offer a potential for blocking cytokine storm-related immunopathology.[40,43,44]

Furthermore, literature review on the most enriched compounds in the compound–target–pathway network indicated their strong correlation to the identified molecular targets in in vitro and in vivo studies. For instance, quercetin was shown to play a modulating, biphasic and regulatory action on inflammation and immunity as it could inhibit lipopolysaccharide (LPS)-induced mRNA levels of TNF-α and interleukin (IL)-1α. Quercetin effectively mitigated rhinovirus-induced progression of lung disease in a mouse model of COPD through prevention of accumulation of neutrophils, macrophages and T cells in the lungs.[45] Quercetin inhibited serum necrosis factor α, interleukin 1β, and interleukin 6, and nitric oxide (NO) in LPS-induced acute lung injury,[46] and suppressed LPS-induced lung inflammation through a heme oxygenase-1-dependent pathway.[47] Similar suggested mechanisms were shown for the flavonoids kaempferol, isorhamnetin, luteolin and apigenin.

Ursolic acid showed inhibitory effect on pulmonary tissues damage in mice through suppression of inflammatory cytokine and oxidative enzymes.[48] It reduced inflammatory factors TNF-α, IL-6 and IL-1β secretion in macrophages in response to LPS stimulation.[49] β glycyrrhetinic acid, isoliquiritigenin, and ursolic acid inhibited the gene expressions of TNF-α, and iNOS, partly through inhibiting NF-κB expression and attenuating NF-κB nuclear translocation.[50]

In addition, glycyrrhizin was revealed to have anti-inflammatory and protective effects on LPS-induced ALI in mice through inhibition of pro-inflammatory cytokines playing a key role in the initial phase of inflammatory response of pulmonary inflammation, which suggests that inhibition of the TLR-4/NF-κB signal pathway would be a possible mechanism underlying its action.[51]

For the identification of the related pathways and GO functional analyses, the identified 26 targets were forwarded to DAVID online database. The obtained results were utilized for the construction of compound–target–pathway network (Table 3, Fig. 8). It could be observed that enriched pathways of the molecular targets were cytokine–cytokine receptor interaction, TNF signaling pathway, NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, NF-kappa B (Nuclear Factor kappa-light-chain-enhancer of activated B cells) signaling pathway, pertussis, JAK-STAT3 signaling pathway, influenza A, tuberculosis, and HTLV-I infection.

KEGG pathway analysis of the identified pathways in the network

# term ID	Term description	Observed gene count	False discovery rate	Matching proteins in network
hsa04668	TNF signaling pathway	11	1.30 × 10−17	AKT1, CASP3, CASP8, FOS, IL1B, JUN, MAPK8, MMP9, PTGS2, TNF, VCAM1
hsa04657	IL-17 signaling pathway	10	2.27 × 10−16	CASP3, CASP8, FOS, HSP90AA1, IL1B, JUN, MAPK8, MMP9, PTGS2, TNF
hsa05161	Hepatitis B	10	1.08 × 10−14	AKT1, CASP3, CASP8, CASP9, FOS, JUN, MAPK8, MMP9, STAT1, TNF
hsa05152	Tuberculosis	10	6.08 × 10−14	AKT1, CASP3, CASP8, CASP9, IL10, IL1B, MAPK8, NOS2, STAT1, TNF
hsa04010	MAPK signaling pathway	9	2.33 × 10−10	AKT1, CASP3, FOS, IL1B, JUN, KDR, MAPK8, TNF, VEGFA
hsa05133	Pertussis	8	3.76 × 10−13	CASP3, FOS, IL10, IL1B, JUN, MAPK8, NOS2, TNF
hsa04620	Toll-like receptor signaling pathway	8	3.94 × 10−12	AKT1, CASP8, FOS, IL1B, JUN, MAPK8, STAT1, TNF
hsa04659	Th17 cell differentiation	8	3.94 × 10−12	FOS, HSP90AA1, IL1B, IL2, JUN, MAPK8, SMAD2, STAT1
hsa04210	Apoptosis	8	2.75 × 10−11	AKT1, CASP3, CASP8, CASP9, FOS, JUN, MAPK8, TNF
hsa05168	Herpes simplex infection	8	2.33 × 10−10	CASP3, CASP8, FOS, IL1B, JUN, MAPK8, STAT1, TNF
hsa05166	HTLV-I infection	8	2.31 × 10−9	AKT1, FOS, IL2, JUN, MAPK8, SMAD2, TNF, VCAM1
hsa04151	PI3K-Akt signaling pathway	8	2.24 × 10−8	AKT1, CASP9, HSP90AA1, IL2, KDR, NOS3, PIK3CG, VEGFA
hsa05321	Inflammatory bowel disease (IBD)	7	1.14 × 10−11	IL10, IL1B, IL2, JUN, SMAD2, STAT1, TNF
hsa04621	NOD-like receptor signaling pathway	7	4.91 × 10−9	CASP8, HSP90AA1, IL1B, JUN, MAPK8, STAT1, TNF
hsa05164	Influenza A	7	5.12 × 10−9	AKT1, CASP9, IL1B, JUN, MAPK8, STAT1, TNF
hsa05165	Human papillomavirus infection	7	2.72 × 10−7	AKT1, CASP3, CASP8, PTGS2, STAT1, TNF, VEGFA
hsa04370	VEGF signaling pathway	6	7.66 × 10−10	AKT1, CASP9, KDR, NOS3, PTGS2, VEGFA
hsa04660	T cell receptor signaling pathway	6	1.09 × 10−8	AKT1, FOS, IL10, IL2, JUN, TNF
hsa04060	Cytokine–cytokine receptor interaction	6	1.85 × 10−6	IL10, IL1B, IL2, KDR, TNF, VEGFA
hsa04658	Th1 and Th2 cell differentiation	5	2.91 × 10−7	FOS, IL2, JUN, MAPK8, STAT1
hsa05144	Jak-STAT signaling pathway	4	1.28 × 10−6	AKT1, IL10, IL2, STAT1
hsa04064	NF-kappa B signaling pathway	4	1.41 × 10−5	IL1B, PTGS2, TNF, VCAM1
hsa04066	HIF-1 signaling pathway	4	1.66 × 10−5	AKT1, NOS2, NOS3, VEGFA
hsa04068	FoxO signaling pathway	4	4.41 × 10−5	AKT1, IL10, MAPK8, SMAD2
hsa05160	Hepatitis C	4	4.41 × 10−5	AKT1, MAPK8, STAT1, TNF
hsa04630	Malaria	4	8.87 × 10−5	IL10, IL1B, TNF, VCAM1
hsa05169	Epstein–Barr virus infection	4	0.00018	AKT1, IL10, JUN, MAPK8
hsa04014	Ras signaling pathway	4	0.00029	AKT1, KDR, MAPK8, VEGFA
hsa04662	B cell receptor signaling pathway	3	0.00018	AKT1, FOS, JUN
hsa04920	Adipocytokine signaling pathway	3	0.00018	AKT1, MAPK8, TNF
hsa01521	EGFR tyrosine kinase inhibitor resistance	3	0.00024	AKT1, KDR, VEGFA
hsa04611	Platelet activation	3	0.00078	AKT1, NOS3, PIK3CG
hsa05310	Asthma	2	0.001	IL10, TNF
hsa04672	Intestinal immune network for IgA production	2	0.0022	IL10, IL2
hsa00590	Arachidonic acid metabolism	2	0.0038	ALOX5, PTGS2
hsa04750	Inflammatory mediator regulation of TRP channels	2	0.0078	IL1B, MAPK8
hsa05322	Systemic lupus erythematosus	2	0.0081	IL10, TNF
hsa04670	Leukocyte transendothelial migration	2	0.0111	MMP9, VCAM1
hsa04062	Chemokine signaling pathway	2	0.0236	AKT1, STAT1

Fig. 8

(A) GO enrichment analysis identified targets. Biological processes are colored orange, cellular components are blue and molecular functions are green. (B) BBID (blue), BIOCARTA (green), INTERPRO (grey) and KEGG (orange) pathways analysis involved in inflammation and immunity. The order of importance was ranked by −log10(P-value) with bar chart. The number of target stick into each term with line chart.

The identified pathways are all closely related to inflammatory, innate and adaptive immune responses. TNF signaling pathway plays an important role in the regulation of immune cells. It is capable of inhibiting viral replication via IL1 and IL6 producing cells. Meanwhile, Toll-like receptors (TLRs) and Nod-like receptors (NLRs) are important mediators of immune recognition. They initiate innate immune responses as well as activation of adaptive immune responses through identification of pathogen-associated molecular patterns. NOD (nucleotide-binding oligomerization domain)-like receptor signaling pathway is responsible for detecting specific pathogen-associated molecules or host-derived damage signals in the cytosol and initiating the innate immune response. NF-kB signaling pathway is involved in stress responses and is activated by pro-inflammatory cytokines like IL-1b and TNF-a. JAK-STAT3 signaling pathway (Janus kinase (JAK)-signal transducer and activator of transcription) plays critical roles in organizing the immune system, especially cytokine receptors and is involved in the modulation of the polarization of T helper cells.

Most researches have shown that SARS-CoV-2 infection can lead to uncontrolled inflammatory innate responses and impaired adaptive immune responses resulting in harmful tissue damage. Drugs directed against key inflammatory cytokines or other aspects of the innate immune response are considered a potential class of adjunctive therapies for SARS-CoV-2 infection.[52,53] Corona virus infection of monocytes and macrophages results in their activation and secretion of IL-6 and other inflammatory cytokines which downstream signal transduction is mediated by JAKs and STAT3.[41,54-56]

Furthermore, results of GO enrichment analysis (Fig. 8) revealed that these targets are involved in a variety of biological processes, including inflammatory response, positive regulation of fever generation, positive regulation of NF-kappaB import into nucleus, positive regulation of interleukin-6 production, positive regulation of interferon-gamma production, prostaglandin biosynthetic process, leukotriene metabolic process, positive regulation of T cell and B cell proliferation in the cellular components extracellular region, secretory granule, organelle membrane and membrane raft.

The 5 enriched molecular functions of the target proteins were mainly associated with cytokine activity, growth factor activity, prostaglandin-endoperoxide synthase activity and heme binding.

The dissemination of the 1073 compound–target interactions on the plants was then inspected (Table S1†) and the plants were ordered in terms of their combined score of compound–target interactions and a merged plant-constituent-target–pathway network was constructed for the top interacting plants (Fig. 7A). The plants with the highest percentage of C–T interactions were Glycyrrhiza glabra (liquorice), Hibiscus sabdariffa, Cichorium intybus, Chrysanthemum coronarium, Nigella sativa, Anastatica hierochuntica, Euphorbia species, Psidium guajava and Epilobium hirsutum.

The plant Glycyrrhiza glabra (liquorice) has a long history of use for the treatment of several inflammatory disease and its effect has been verified through different in vitro and in vivo models. In one study, Glycyrrhiza glabra suppressed the ability of LPS to induce mRNA and protein expressions of tumor necrosis factor-α (TNF-α), interleukin 1β (IL-1β), and IL-6 productions.[57] In another study, Glycyrrhiza glabra displayed preventive and therapeutic effects on pulmonary fibrosis.[58] The flavones of liquorice markedly decreased pro-inflammatory cytokines expression levels through modulation of NF-kB/MAPK pathway.[59]

Meanwhile, Hibiscus sabdariffa was revealed to possess potent anti-inflammatory activity as it maintained the ratio of IL-1β/IL-1ra in the plasma and hippocampus of Wistar rats model.[60] In another study, Hibiscus sabdariffa and its isolate delphinidin in a cell model reduced the levels of inflammatory mediators including iNOS, IL-6, and TNF-α induced by LPS and downregulated NF-κB pathway. In animal model, they reduced the production of IL-6, and TNF-α and attenuated mouse paw edema induced by LPS.[61]

Cichorium intybus L. (chicory) is commonly used for the treatment of inflammatory conditions. Chicory roots demonstrated significant dose-dependent decrease in paw edema in carrageenan-induced paw edema method. Chicory roots diminished the serum TNF-α, IL-6, and IL-1 levels[62] and its anti-inflammatory action was shown to be due to a direct modulation of cytokine expression.[63] The plant Chrysanthemum coronarium was shown to decrease TNF-α synthesis in the cells of monocyte origin activated with LPS, and decreased cytokine levels in vitro.[64]

The plant Nigella sativa has a long history of use as an immunomodulatory and anti-inflammatory herb. In an in vitro study, it differently modulated the pro-inflammatory cytokines including Interleukin (IL)-1 alpha, IL-1 beta, IL-6.[65]Nigella sativa also is known to possess the potential to modulate adaptive immunity and may play a critical role in modulating the balance of Th1/Th2 lymphocytes, leading to altered Th1/Th2 cytokine profiles.[66]

Conclusion

Many of the medicinal plants and their constituents have a potential for use in the mitigation of the new SARS-CoV-2 infection. Herein, a database comprised of more than 16 500 compounds was screened against the three viral targets 3CLpro, PLpro and RdRp, and several constituents identified may inhibit SARS-CoV-2 activity through inhibition of virus replication. A network pharmacology analysis was performed for all the plants constituents and revealed that several compounds possessed multi-targets and enriched pathways of the molecular targets including cytokine–cytokine receptor interaction, TNF signaling pathway, Toll-like receptor signaling pathway, NF-kappa B signaling pathway, and JAK-STAT3 signaling pathway. These results suggest a potential role of medicinal plants in the management of the current SARS-CoV-2 infection. A special emphasis is paid to the plant Glycyrrhiza glabra (liquorice) which has a long history of use in treatment of several viral infections. In fact, randomized controlled trials confirmed that the Glycyrrhiza glabra demonstrated a reduction of mortality and viral activity in SARS-CoV-2 related coronavirus.[67] Meanwhile, the plants Hibiscus sabdariffa and Cichorium intybus are rich sources of caffeic acid derivatives which have been widely investigated concerning their antiviral potential.[68-70]

The performed network pharmacology analysis revealed the synergistic nature of the compounds within each medicinal plant as shown in Fig. 6 and 7. Combination of the most enriched plants in the created network i.e. Glycyrrhiza glabra (liquorice), Hibiscus sabdariffa, Cichorium intybus, could be chosen for mitigation of SARS-CoV-2.

Further molecular dynamic simulation studies, in vitro and in vivo analyses are required to confirm the results of our study, however, the insights provided by the present study may substantiate valuable exploration and development of anti-SARS-CoV-2 therapeutic agents from natural origin.

Conflicts of interest

None.