Analysis of mechanisms of Shenhuang Granule in treating severe COVID-19 based on network pharmacology and molecular docking.

OBJECTIVE: Severe cases of coronavirus disease 2019 (COVID-19) are expected to have a worse prognosis than mild cases. Shenhuang Granule (SHG) has been shown to be a safe and effective treatment for severe COVID-19 in a previous randomized clinical trial, but the active chemical constituents and underlying mechanisms of action remain unknown. The goal of this study is to explore the chemical basis and mechanisms of SHG in the treatment of severe COVID-19, using network pharmacology.
METHODS: Ultraperformance liquid chromatography/quadrupole time-of-flight mass spectrometry was employed to screen chemical constituents of SHG. Putative therapeutic targets were predicted by searching traditional Chinese medicine system pharmacology database and analysis platform, SwissTargetPrediction, and Gene Expression Omnibus (GEO) databases. The target protein-protein interaction network and enrichment analysis were performed to investigate the hub genes and presumptive mechanisms. Molecular docking and molecular dynamics simulations were used to verify the stability and interaction between the key chemical constituents of SHG and COVID-19 protein targets.
RESULTS: Forty-five chemical constituents of SHG were identified along with 131 corresponding therapeutic targets, including hub genes such as HSP90AA1, MMP9, CXCL8, PTGS2, IFNG, DNMT1, TYMS, MDM2, HDAC3 and ABCB1. Functional enrichment analysis indicated that SHG mainly acted on the neuroactive ligand-receptor interaction, calcium signaling pathway and cAMP signaling pathway. Molecular docking showed that the key constituents had a good affinity with the severe acute respiratory syndrome coronavirus 2 protein targets. Molecular dynamics simulations indicated that ginsenoside Rg4 formed a stable protein-ligand complex with helicase.
CONCLUSION: Multiple components of SHG regulated multiple targets to inhibit virus invasion and cytokine storm through several signaling pathways; this provides a scientific basis for clinical applications and further experiments.

Introduction

The pandemic of coronavirus disease 2019 (COVID-19) has caused worldwide concern. As of November 2021, there have been over 248 million confirmed cases and more than 5 million deaths worldwide [1]. An early cohort study of patients with COVID-19 showed that 14% had severe disease, 5% had critical disease, and the mortality of the critically ill cases was 49 % [2]. Severe COVID-19 cases are expected to have a worse prognosis than milder cases. A systematic review of global COVID-19 infections showed that 67.7% of patients admitted to intensive care units required invasive mechanical ventilation, 65.9 % required vasopressor support, 16.9% required renal replacement therapy, and 6.4% required extracorporeal membrane oxygenation [3]. Therapeutic considerations for severe and critical COVID-19 include supportive care, corticosteroids, remdesivir and tocilizumab [4]. Although remarkable advances have been made in reducing mortality and enhancing recovery from symptoms, no drugs have been developed to treat the infection of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2).

Traditional Chinese medicine (TCM) has gained important experience in managing the COVID-19 epidemic in China. According to the Diagnosis and Treatment Protocol for Novel Coronavirus Pneumonia (8th edition), Qingfei Paidu Decoction, Xuanfei Paidu Formula, Huashi Baidu Formula, Jinhua Qinggan Granule, Lianhua Qingwen Capsule and Xuebijing Injection, referred to as “three Chinese patent medicines and three decoctions,” are recommended for treating the disease [5]. Among these, Xuebijing Injection and Huashi Baidu Formula are used mainly for severe COVID-19 [6], [7], while Qingfei Paidu Decoction can be used for patients at any clinical stage [8]. In the early stage of outbreaks, our group developed the Shenhuang Granule (SHG) to treat severe COVID-19. We also conducted a multicenter, prospective and randomized clinical trial to confirm the efficacy and safety of SHG in patients with severe COVID-19 [9]. The result reflected an increased recovery rate and a decreased mortality in severe and critical cases. Further, SHG was able to reduce the rate of disease progression from the severe category to the critical category.

SHG is composed of Radix Ginseng, Radix Aconiti Lateralis Preparata, Rhei Radix Et Rhizoma, Caulis Sargentodoxae, Herba Taraxaci and Hirudo. As a complex herbal mixture containing multiple chemical constituents that act on numerous targets and pathways, SHG may be effective against COVID-19. However, the active chemical constituents, putative targets and biological mechanisms of SHG remain unknown and need further investigation. Network pharmacology is a new subject based on systems biology, computational biology and systems pharmacology. Network pharmacology explores the pharmacological characteristics and mechanisms of TCM prescriptions in the treatment of diseases from a systematic and comprehensive point of view by constructing a “disease-function-target-constituent-herb” interaction network, analyzing the topological parameters, and calculating the relevance between the specific signal nodes in the network [10]. Then, the key chemical constituents in the network and the core targets of the disease are selected for molecular docking simulation, in order to obtain the potential effective chemical constituents and targets for further experimental verification. Recent network pharmacology studies have preliminarily illustrated the interaction mechanisms between small molecules from Chinese medicines and COVID-19 disease targets. For example, quercetin and baicalein, present in the Huashi Baidu Formula, were found to have a strong affinity for angiotensin-converting enzyme 2 (ACE2) and SARS-CoV-2 3-chymotrypsin-like protease protein [11]. Baicalein, luteolin and quercetin in the Xuebijing Injection may inhibit the excessive inflammatory response associated with COVID-19 by regulating AKT1 expression [12]. Four compounds (baicalin, glycyrrhizic acid, hesperidin and hyperoside) in the Qingfei Paidu Decoction have interactions with different SARS-CoV-2 targets [13]. These results provide a new rationale for Chinese medicine to treat COVID-19.

In this study, we identified the main chemical constituents of SHG by ultraperformance liquid chromatography/quadrupole time-of-flight mass spectrometry (UPLC-Q-TOF-MS). Putative therapeutic targets were predicted from several databases. Then, network pharmacology methods, molecular docking, and molecular dynamics simulations were integrated to investigate the mechanisms of SHG on severe COVID-19 and verify the stability and interaction between the key chemical constituents and COVID-19 protein targets.

Materials and methods

Materials and chemicals

Liquid chromatography-mass spectrometry (LC-MS) grade methanol and acetonitrile were provided by Merck KGaA (Shanghai, China). The same grade methanoic acid (Y9330090) was purchased from Thermo Fisher Scientific Co., ltd (Shanghai, China). The deionized water (20210402C) was provided by Guangzhou Watsons Food & Beverage Co., Ltd (Guangzhou, China).

Preparation of SHG extraction

This study investigated the Chinese patent medicine SHG, which is composed of Radix Ginseng (Renshen) 50 g, Rhei Radix Et Rhizoma (Dahuang) 40 g, Radix Aconiti Lateralis Preparata (Fuzi) 50 g, Caulis Sargentodoxae (Daxueteng) 30 g, Herba Taraxaci (Pugongying) 30 g, and Hirudo (Shuizhi) 6 g. SHG was produced by Beijing Tcmages Pharmaceutical Co. ltd (Beijing, China). Quality assurance test results for each drug were consistent with the required quality standards. The product was approved by the National Medical Product Administration (China; Approval number: Jing 20180032).

UPLC-Q-TOF-MS identification of major SHG constituents

A 0.5 g sample of SHG was added to 25 mL of 80 % methanol in a 50 mL corked conical flask. The SHG was solubilized with ultrasound (300 W, 40 kHz) for 30 min. Then the SHG solution was centrifuged at 13,000 × g for 5 min to obtain supernatant. Analysis was performed on an Agilent 1290 UPLC system (Agilent, USA) in tandem with an AB Sciex Triple TOF® 4600 mass spectrometer (AB SCIEX, USA). Experimental parameters were set according to Hao et al [14]. UPLC parameters were as follows: a Waters Acquity UPLC HSS T3 column (2.1 mm × 100 mm, 1.8 μm), with a flow rate of 0.3 mL/min, column temperature 30 ℃, and detection wavelength from 190 − 400 nm. Mobile phase A was 0.1% methanoic acid in water, and phase B was 0.1 % methanoic acid in acetonitrile. The gradient elution protocol was as follows: 0–14 min, 5%−20% B; 14–39 min, 20%−60% B; 39–46 min, 60%−95% B; 46.1–48 min, 95%−5% B. Mass spectrometry (MS) conditions were as follows: time-of-flight mass range, 50–1700; ion source gas 1, 2.394 kPa; ion source gas 2, 2.394 kPa; curtain gas, 1.676 kPa; ion spray voltage floating, –4500/5000 V; ion source temperature, 500 °C; declustering potential, 100 V; MS/MS mass range, 50–1250; collision energy, ± 40 eV. MS data were collected using an Analyst TF 1.7.1 (AB SCIEX, Foster City, CA, USA) and processed using Peakview 1.2 (AB SCIEX, USA). The main peaks in the fingerprints were identified and authenticated using the Natural Products HR-MS/MS Spectral Library 1.0 software (Shanghai Standard Technology Co., ltd., Shanghai, China), which contains over 1000 standard substances. Chemical constituents that were not found in the database were identified using literature reports and the fragmentation patterns in MS.

Screening of putative targets in SHG

The traditional Chinese medicine system pharmacology database and analysis platform (TCMSP, https://tcmsp-e.com/) [15], and the SwissTargetPrediction database (https://www.swisstargetprediction.ch/) [16] were used to identify the putative targets of the main constituents of SHG. The protein target corresponding to each constituent was translated into a gene name using the UniProt database (https://www.uniprot.org/).

Screening of severe COVID-19 disease targets

The gene expression array datasets of severe COVID-19 were searched using the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) [17]. Microarray data for differentially expressed long noncoding RNAs (lncRNAs) in the peripheral blood mononuclear cells between the severe COVID-19 group and the healthy group were obtained from GEO Series GSE164805 (PMID: 33679778) [18]. The GEO2R online tool was used to detect the genes that are differentially expressed between the two groups, and genes with an adjusted P < 0.05 and |log2(fold change)| > 2, were considered significantly differentially expressed and potential drug targets in severe COVID-19. Finally, the volcano map of significantly differentially expressed genes was visualized using the NetworkAnalyst platform (https://www.networkanalyst.ca/NetworkAnalyst/home.xhtml) [19].

Construction of SHG herb-constituent-target network

The therapeutic targets of SHG in the treatment of severe COVID-19 were obtained by intersecting the putative targets of the chemical constituents with the disease targets. A Venn diagram was then prepared using the Venn diagram webtool (https://bioinformatics.psb.ugent.be/cgi-bin/liste/Venn/). The SHG herb-constituent-target network was constructed to illustrate the relationship between the herb chemical constituents and the severe COVID-19 disease targets. Cytoscape 3.8.2 software [20] was used to construct the herb-constituent-target network, and the network degree parameter was calculated. Deeper colors and larger node areas represent greater connectivity among targets and chemical constituents.

Construction of the protein–protein interaction network and identification of hub genes

The protein–protein interaction (PPI) network of the therapeutic targets was obtained using the STRING database (https://string-db.org/) [21]. The therapeutic targets were input into STRING to remove the poorly connected targets and the threshold of an interaction score of > 0.4 was used. The PPI network was visualized using Cytoscape 3.8.2. The CytoNca plugin was used to calculate the degree of each protein node and screen the hub genes. Three important parameters that represented the topological importance of the nodes in the network were selected, including betweenness centrality, closeness centrality and degree centrality. We selected the top ten genes as the hub genes. The relationship between hub genes and chemical constituents of SHG is shown through the Sankey diagram.

Construction of target-organ network

Considering that COVID-19 is a systemic disease, multiple tissues and organs may suffer different degrees of damage over the course of severe COVID-19. We used the BioGPS database (https://biogps.org) [22] to analyze the tissue-specific expression and the distribution of the 131 therapeutic targets. The screening criteria were as follows: the tissue-specific expression level was > 10 times the median; if not, the top 3 expression levels were extracted. The constituent target-organ location network was visualized using Cytoscape 3.8.2. The number of adjacent nodes in the network was calculated, and the size of the nodes in the network was determined according to the number of adjacent nodes.

Gene Ontology and Kyoto Encyclopedia of genes and Genomes enrichment analysis

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were applied to the annotated target proteins at the functional level. The ClusterProfiler package [23] of R 3.6.3 was used to analyze the GO and KEGG enrichment. An adjusted P < 0.05 and a q value < 0.2 were considered as a threshold for restriction. The top 20 GO pathways in the biological process (BP), cellular composition (CC) and molecular function (MF) as well as the top 20 KEGG pathways were retained.

Physicochemical property and ADMET evaluation of the main chemical constituents in SHG

Based on the MolSoft website (https://molsoft.com/mprop/), the drug-likeness scores of the main chemical constituents in SHG were calculated by querying simplified molecular input line entry system. Based on the ADMETlab 2.0 server (https://admetmesh.scbdd.com/) [24], the physicochemical properties of each compound were predicted, including molecular weight, number of hydrogen bond acceptors, number of hydrogen bond donors, LogP, LogS and number of stereo centers. In addition, the ADMET parameters were predicted, namely absorption, distribution, metabolism, excretion and toxicity.

Molecular docking

The key chemical constituents with more than the median degree, betweenness centrality and closeness centrality were selected from the herb-constituent-target network, and the SDF files of these constituents were downloaded from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) and then uploaded to the COVID-19 Docking Server website (https://ncov.schanglab.org.cn/index.php) [25] for molecular docking. The COVID-19 Docking Server website serves as a platform for docking small molecules, peptides, or antibodies to COVID-19 protein targets and provides a prediction of target-ligand interactions. The protein targets include 1 structural protein (nucleocapsid protein [N protein]) and 8 nonstructural proteins (main protease [Mpro], papain-like protease [PLpro], nonstructural protein 3 [Nsp3], RNA-dependent RNA polymerase [RdRp], helicase, Nsp14, Nsp15 and Nsp16) that are considered ideal targets for developing anti-viral drugs. The Nsp3 has two binding sites, referred to as adenosine monophosphate (AMP) and 2-(N-morpholino)-ethanesulfonic acid sites. The RdRp has two binding sites referred to as RdRp without RNA and RdRp with RNA. The helicase has two binding sites, called adenosine diphosphate and nucleic acid binding sites. The Nsp14 has two binding sites, called N-terminal exoribonuclease and guanine-N7 methyl transferase sites. The Nsp16 has three binding sites, called 7-methyl-GpppA (GTA), S-adenosylmethionine (SAM) and 7-methyl-guanosine-5′-triphosphate sites. The corresponding Protein Data Bank (PDB) codes were 6LU7 (Mpro), 6WUU (PLpro), 6W6Y (Nsp3), 7BV2 (RdRp), 6JYT (helicase), 5C8S (Nsp14), 6WLC (Nsp15), 6WVN (Nsp16) and 4KXJ (N protein). PyMol (Version 2.5.0, Schrödinger, Germany) and LigPlot (Version 2.2.5, EMBL’s European Bioinformatics Institute, United Kingdom) software were used to examine the binding mode and interactions between the selected chemical constituents and the COVID-19 protein targets.

Molecular dynamics simulation

The molecular dynamics simulation was used to verify the stability and dynamic interactions between the proteins and ligands after molecular docking. The molecular dynamics simulation was performed using the GROMACS 2018.8 software (https://www.gromacs.org/), running on the CentOS Linux 7 operating system (Red Hat, USA). The Charmm36 all-atom force field [26] was used to process the protein system and the CGenFF server (https://cgenff.umaryland.edu) [27] was used to generate the ligand topology. The system was solvated using the TIP3P water model in a dodecahedron box, and appropriate numbers of Na+/Cl− ions were added to neutralize charges. The system was energy-minimized using the steepest descent method of 50,000 steps and < 10.0 kJ/mol force and subsequently equilibrated using 100 ps of canonical and isothermal-isobaric simulations at the temperature of 300 K. Finally, the molecular dynamics simulation was carried out using a leap-frog algorithm for 10 ns with a time step of 2 fs. The resulting molecular dynamics trajectories were analyzed for several parameters including root-mean-square deviation and root-mean-square fluctuation. The Gibbs free energy was calculated according to the values of root-mean-square deviation and radius of gyration (Rg). Xmgrace 5.1.22 (https://plasma-gate.weizmann.ac.il/Grace/) and Origin 8.5.1 (OriginLab, USA) were used for visualization.

Results

Information on chemical constituents and putative targets of SHG

A total of 45 main chemical constituents of SHG were detected by UPLC-Q-TOF-MS (Table 1 ). The base peak chromatograms of SHG under both positive- and negative-ion modes are shown in Fig. 1 . Among these chemical constituents, we could not determine the specific molecular structure for peaks 41, 43 and 18. Further, components numbered 24 and 36 were not found on the PubChem website and were drawn from the SwissTargetPrediction database. Searches of the TCMSP and SwissTargetPrediction databases returned 829 targets for the 42 chemical constituents, after removing duplicate values.

Table 1

Main chemical constituents identified from Shenhuang Granule using ultraperformance liquid chromatography/quadrupole time-of-flight mass spectrometry.

No.	Measured mass (m/z)	Error (‰)	Molecular formula	Molecule name	Herb	PubChem CID
1	331.0671	0.0030	C13H16O10	O-Galloylglucose	Rhei Radix Et Rhizoma	21120349
2	169.0142	0.0070	C7H6O5	Gallic acid	Rhei Radix Et Rhizoma	370
3	164.0717	0.0049	C9H11NO2	Phenprobamate	Hirudo	6140
4	353.0878	0.0010	C16H18O9	Neochlorogenic acid	Herba Taraxaci/Caulis Sargentodoxae	5280633
5	486.2698	–0.0051	C24H39NO9	Mesaconine	Radix Aconiti Lateralis Preparata	76189547
6	299.1136	0.0026	C14H20O7	Salidroside	Caulis Sargentodoxae	159278
7	289.0718	0.0036	C15H14O6	Catechin	Herba Taraxaci/Caulis Sargentodoxae/Rhei Radix Et Rhizoma	1203
8	353.0878	0.0040	C16H18O9	Chlorogenic acid	Herba Taraxaci/Caulis Sargentodoxae	1794427
9	353.0878	0.0020	C16H18O9	Cryptochlorogenic acid	Herba Taraxaci/Caulis Sargentodoxae	9798666
10	454.2799	–0.0036	C24H39NO7	Fuziline	Radix Aconiti Lateralis Preparata	14163819
11	787.2666	0.0020	C34H46O18	Eleutheroside E	Caulis Sargentodoxae	71312557
12	473.0725	0.0010	C22H18O12	Chicoric acid	Herba Taraxaci	5281764
13	431.0984	0.0040	C21H20O10	Aloe-emodin-8-O-ß-D-glucopyranoside	Rhei Radix Et Rhizoma	5317644
14	445.0776	0.0030	C21H18O11	Rhein-8-O-ß-D-glucopyranoside	Rhei Radix Et Rhizoma	348570906
15	590.2960	–0.0030	C31H43NO10	Benzoylmesaconine	Radix Aconiti Lateralis Preparata	24832659
16	431.0984	0.0024	C21H20O10	Emodin-1-O-glucoside	Rhei Radix Et Rhizoma	5319333
17	604.3116	–0.0030	C32H45NO10	Benzoylaconine	Radix Aconiti Lateralis Preparata	20055771
18*	459.0933	0.0022	C22H20O11	3(or 8)-(ß-D-Glucopyranosyloxy)-9,10-dihydro-8(or 3)-hydroxy-1-methyl-9,10-dioxo-2-anthracenecarboxylic acid	Rhei Radix Et Rhizoma	-
19	574.3011	–0.0029	C31H43NO9	Benzoylhypaconine	Radix Aconiti Lateralis Preparata	78358526
20	579.2083	0.0030	C27H34O11	Arctiin	Herba Taraxaci	100528
21	845.4904	0.0030	C42H72O14	Ginsenoside Rg1	Radix Ginseng	441923
22	991.5483	0.0030	C48H82O18	Ginsenoside Re	Radix Ginseng	441921
23	616.3116	–0.0013	C33H45NO10	Hypaconitine	Radix Aconiti Lateralis Preparata	441737
24	588.3167	–0.0030	C32H45NO9	3-Deoxyhokbusine A	Radix Aconiti Lateralis Preparata	-
25	415.1035	0.0023	C21H20O9	Chrysophanol-8-O-ß-D-glucopyranoside	Rhei Radix Et Rhizoma	442731
26	431.0984	0.0024	C21H20O10	Aloe-emodin-3-(hydroxymethyl)-O-ß-D-glucopyranoside	Rhei Radix Et Rhizoma	147295
27	431.0984	0.0031	C21H20O10	Emodin-8-O-glucoside	Rhei Radix Et Rhizoma	5319333
28	445.1140	0.0087	C22H22O10	Physcion 8-O-ß-D-monoglucoside	Rhei Radix Et Rhizoma	168938
29	845.4904	0.0030	C42H72O14	Ginsenoside Rf	Radix Ginseng	441922
30	269.0455	0.0028	C15H10O5	Aloe-emodine	Rhei Radix Et Rhizoma	10207
31	1153.6011	0.0050	C54H92O23	Ginsenoside Rb1	Radix Ginseng	9898279
32	955.4908	0.0050	C48H76O19	Ginsenoside Ro	Radix Ginseng	11815492
33	1123.5906	0.0050	C53H90O22	Ginsenoside Rb2	Radix Ginseng	6917976
34	283.0248	0.0007	C15H8O6	Rhein	Rhei Radix Et Rhizoma	10168
35	991.5483	0.0030	C48H82O18	Ginsenoside Rd	Radix Ginseng	11679800
36	589.1352	0.0025	C31H26O12	2-Anthracenecarboxylic acid, 9,10-dihydro-8-hydroxy-1-methyl-9,10-dioxo-3-[[6-O-[(2E)-1-oxo-3-phenyl-2-propen-1-yl]-ß-D-glucopyranosyl]oxy]-	Rhei Radix Et Rhizoma	-
37	811.4849	0.0014	C42H70O12	Ginsenoside Rg4	Radix Ginseng	102004835
38	665.4270	0.0040	C36H60O8	Ginsenoside Rh4	Radix Ginseng	21599928
39	269.0455	0.0028	C15H10O5	Emodin	Rhei Radix Et Rhizoma	3220
40	793.4380	0.0050	C42H66O14	Chikusetsusaponin Iva	Radix Ginseng	13909684
41*	829.4955	0.0030	C42H72O13	Ginsenoside Rg3 isomer	Radix Ginseng	-
42	253.0506	0.0070	C15H10O4	Chrysophanol	Rhei Radix Et Rhizoma	10208
43*	829.4955	0.0020	C42H72O13	Ginsenoside Rg3 isomer	Radix Ginseng	-
44	811.4849	0.0020	C42H70O12	Ginsenoside Rk1	Radix Ginseng	11499198
45	811.4849	0.0020	C42H70O12	Ginsenoside Rg5	Radix Ginseng	11550001

*: Chemical constituents without a specific molecular structure. Numbers in the No. column are the same as the peak numbers in Fig. 1.

Fig. 1

The base peak chromatograms of Shenhuang Granule. A: positive-ion mode; B: negative-ion mode. The details of the constituents represented by the peak numbers (clustered peaks marked with black arrows) are shownin Table 1.

Potential target genes of severe COVID-19

The analysis of the gene chip in the GEO database (GSE164805) identified 6160 differentially expressed genes related to severe COVID-19, which were used to build a volcano plot (Fig. 2 ). These significantly differentially expressed genes are marked with red dots for upregulated genes and blue dots for downregulated genes.

Fig. 2

Volcano plot of differentially expressed genes related to severe coronavirus disease 2019. The red dots indicate upregulated genes and the blue dots indicate downregulated genes. |log2(fold change)| > 2, adjusted P < 0.05.

Herb-constituent-target network construction

The constituent targets of SHG were matched with the disease targets of severe COVID-19, resulting in the selection of 131 therapeutic targets (Fig. 3 ). The herb-constituent-target network of SHG against severe COVID-19 was constructed to illustrate the relationship between the main chemical constituents and the target proteins. The network contained 182 nodes (6 herbs, 45 main chemical constituents, and 131 target proteins) and 761 edges (Fig. 4 ). The red nodes represented the herbs. The peripheral rings were chemical constituents of the herbs. In the inner layer, the blue squares represented target proteins. The dark purple and brown hexagons indicate the constituent that existed in multiple herbs. Supplementary Table S1 shows the main chemical constituents of SHG and their corresponding therapeutic targets.

Fig. 3

Venn diagram of the 131 therapeutic targets. SHG: Shenhuang Granule; COVID-19: coronavirus disease 2019.

Fig. 4

Herb-constituent-target network. DH: Rhei Radix Et Rhizoma; PGY: Herba Taraxaci; DXT: Caulis Sargentodoxae; RS: Radix Ginseng; SZ: Hirudo; FZ: Radix Aconiti Lateralis Preparata.

PPI network construction and hub gene identification

Protein interactions among the therapeutic targets were predicted using the STRING database and visualized using Cytoscape. A total of 108 nodes and 329 edges were involved in the PPI network after removing the unconnected targets (Fig. 5 A). CytoNca was used to calculate the degree of each node. The larger the node area and the redder the color, the more important the target protein is anticipated to be. The top ten genes, heat-shock protein 90 α family class a member 1 (HSP90AA1), matrix metallopeptidase 9 (MMP9), C-X-C motif chemokine ligand 8 (CXCL8), prostaglandin-endoperoxide synthase (PTGS2), interferon γ (IFNG), recombinant DNA methyltransferase 1 (DNMT1), Homo sapiens thymidylate synthetase (TYMS), murine double minute 2 (MDM2), histone deacetylase 3 (HDAC3) and ATP-binding cassette subfamily B member 1 (ABCB1), as evaluated by the degree of interaction in the PPI network, were identified as hub genes (Fig. 5C). The Sankey diagram shows the relationship between hub genes, chemical constituents and herbs (Fig. 5B).

Fig. 5

The protein–protein interaction (PPI) network and the hub genes. A: the PPI network; B: the relationship between hub genes, chemical constituents and herbs; C: the top ten hub genes.

Target-organ location network construction

The tissue-specific expression and distribution of the 131 therapeutic targets were analyzed in BioGPS. PDE11A and TLR9 did not have suitable expression profiles in BioGPS. The expression of the remaining 129 targets in a specific tissue or organ system is shown in Supplementary Table S2. Expression of the target proteins was evaluated in specific tissue or organ systems, including the brain, lymphoid tissues, heart, gastrointestinal tract, endocrine tissues, liver, muscle tissues, whole blood, lung, bronchial epithelial cells, adipocyte, male tissues, placenta, nerve tissues, female tissues, eye, bone marrow and kidney. In addition, the targets were highly related to immunocytes, including 721-B-lymphoblasts, CD34+ cells, CD33+ myeloid cells, CD71+ early erythroid cells, CD105+ endothelial cells, BDCA4+ dendritic cells, CD4+ T cells, CD56+ NK cells, CD8+ T cells, CD14+ monocytes, CD19+ B cells and MOLT-4 cells. Then, the target-organ location network, which contained 162 nodes and 320 edges, was constructed to exhibit the relationship between targets and tissues (Fig. 6 ). The blue squares represent target proteins. The peripheral rings indicate tissue or organ locations. Size of the tissue and organ nodes is relative to the degree of network connectivity.

Fig. 6

Target-organ location network of tissue-specific expression gene in BioGPS.

GO and KEGG enrichment analysis

A total of 1103 BP entries, 44 CC entries and 108 MF entries were obtained through the GO enrichment analysis. The KEGG functional enrichment analysis contained 63 items. We selected the top 20 functional enriched terms in the BP, CC and MF categories to draw a bubble diagram (Fig. 7 A−C; Supplementary Table S3). Based on the clustering analysis, BPs were mainly involved in second-messenger-mediated signaling, response to oxygen levels, response to antibiotic, regulation of synaptic transmission, metabolic process, neutrophil activation, and regulation of membrane potential. The CC domain was mainly enriched in the neuronal cell body, synaptic membrane, external side of plasma membrane, and presynapse. The MF domain was mainly involved in phosphoric ester hydrolase activity, endopeptidase activity, amide and peptide binding, and active transmembrane transporter activity. The top 20 KEGG pathways are shown in Fig. 7D and Supplementary Table S4, including neuroactive ligand-receptor interaction, calcium signaling pathway, cAMP signaling pathway, prostate cancer, purine metabolism, alcoholism, interleukin-17 (IL-17) signaling pathway, thyroid hormone signaling pathway, forkhead box O (FoxO) signaling pathway, dopaminergic synapse, apoptosis, estrogen signaling pathway, fluid shear stress and atherosclerosis, gastric cancer, morphine addiction, small cell lung cancer, toll-like receptor (TLR) signaling pathway, C-type lectin receptor (CLR) signaling pathway, T cell receptor (TCR) signaling pathway, and tumor necrosis factor (TNF) signaling pathway. These pathways involve the immune system, signal transduction, cancers, endocrine and nervous system, cell growth and death, cardiovascular disease, metabolism, and substance dependence. The size of the nodes was decided by the number of associated genes, and the colors from red to blue reflected the P-values from high to low.

Fig. 7

The GO and KEGG enrichment analysis with the 131 targets. A: GO biological process (BP) analysis; B: GO cellular composition (CC) analysis; C: GO molecular function (MF) analysis; D: KEGG pathway analysis. KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology.

Physicochemical property and ADMET evaluation of the main chemical constituents

Physicochemical property statistics of the main chemical constituents in SHG are shown in Supplementary Table S5. Among the 42 chemical constituents in SHG, 33 chemical constituents were predicted to possess a positive drug-likeness score. ADMET parameters, including absorption (Caco-2 permeability, human intestinal absorption, and Madin-Darby canine kidney), distribution (plasma protein binding, volume distribution, blood brain barrier, and fraction unbound in plasma), excretion (clearance, and half-life of a drug), and toxicity (human hepatotoxicity, Ames toxicity, rat oral acute toxicity, skin sensitization, eye irritation, and respiratory toxicity) are shown in Supplementary Table S6.

To further assess the interaction between the chemical constituents and the targets, 10 key chemical constituents according to topological analysis and 15 protein targets from the COVID-19 Docking Server were selected for docking. The topology parameters of the 10 key chemical constituents in the herb-constituent-target network are shown in Table 2 . The docking scores between these key chemical constituents of SHG and SARS-CoV-2 protein targets are shown in Table 3 . Chemical constituents with lower binding energies are usually considered to exhibit a higher binding affinity with the protein targets [28]. The affinity scores of the best six modes, chikusetsusaponin Iva-RdRp with RNA, ginsenoside Rf-RdRp with RNA, ginsenoside Rg4-helicase (nucleic acid binding site), ginsenoside Rh4-PLpro, ginsenoside Rg4-Nsp16 (GTA site), and ginsenoside Rg4-Nsp16 (SAM site), were all greater than –9.8 kcal/mol, which indicates that these chemical constituents had a good affinity with the SARS-CoV-2 protein targets (Fig. 8 ).

Table 2

Topology parameter of the key chemical constituents in the network.

No.	Betweenness centrality	Closeness centrality	Degree	Name
39	0.0883	0.3996	25	Emodin
29	0.0451	0.3961	22	Ginsenoside Rf
37	0.0474	0.3892	22	Ginsenoside Rg4
2	0.0758	0.3892	21	Gallic acid
11	0.0409	0.3909	21	Eleutheroside E
14	0.0626	0.3926	21	Rhein-8-O-β-D-glucopyranoside
40	0.0474	0.3876	21	Chikusetsusaponin IVa
38	0.0449	0.3859	20	Ginsenoside Rh4
42	0.0647	0.3843	20	Chrysophanol
34	0.0411	0.3795	17	Rhein

Table 3

Docking scores between the key chemical constituents of Shenhuang Granule and 2019-nCov protein targets.

Protein target	Affinity score (kcal/mol)
	2	11	14	29	34	37	38	39	40	42
Mpro	–5.5	–7.6	–8.3	–7.7	–7.4	–7.1	–6.4	–7.2	–8.6	–7.2
PLpro	–5.9	–9.7	–9.3	–9.1	–8.4	–9.7	–10.1	–8.3	–8.3	–7.9
Nsp3 (AMP site)	–5.4	–7.6	–7.7	–7.9	–7.6	–7.5	–7.7	–7.6	–8.9	–7.7
Nsp3 (MES site)	–6.3	–7.5	–7.6	–7.6	–8.5	–8.9	–7.6	–7.8	–8.4	–7.7
RdRp with RNA	–7.3	–9.8	–9.2	–10.6	–8.6	–9.6	–9.3	–8.4	–11.9	–8.2
RdRp without RNA	–5.7	–7.7	–7.9	–8.1	–7.6	–7.8	–7.7	–7.1	–8.3	–7.1
Helicase (ADP site)	–5.1	–7.1	–6.9	–6.7	–6.8	–7.7	–7.2	–6.5	–7.8	–6.5
Helicase (NCB site)	–5.5	–7.4	–8.5	–8.4	–8	–10.1	–8.8	–7.3	–9.2	–7.5
Nsp14 (ExoN)	–5.3	–7.3	–7.7	–8.6	–7.5	–9	–7.7	–7	–8.9	–7.3
Nsp14 (N7-MTase)	–6.6	–8.6	–9.5	–8.9	–9.1	–8.5	–8.3	–8.8	–8.3	–9.1
Nsp15 (endoribonuclease)	–5.3	–6.4	–7	–7.5	–7.7	–8.1	–7.2	–7.5	–8.1	–7.6
Nsp16 (GTA site)	–6.1	–8.7	–8	–9.2	–8	–9.9	–8	–8	–9.7	–8.1
Nsp16 (MGP site)	–5.1	–7	–8	–7.9	–7	–7.6	–7.7	–7.1	–8.3	–7.2
Nsp16 (SAM site)	–5.9	–8	–8.1	–8.8	–7.3	–9.9	–8.7	–7.2	–9.6	–7.1
N protein	–5.4	–7.5	–8.8	–7.7	–8.6	–8.7	–8.3	–8.8	–8.7	–8.6

The details of the constituents represented by the numbers are shown in Table 1. ADP: adenosine diphosphate; AMP: adenosine monophosphate; ExoN: N-terminal exoribonuclease; GTA: 7-methyl-GpppA; MES: 2-(N-morpholino)-ethanesulfonic acid; MGP: 7-methyl-guanosine-5′-triphosphate; Mpro: main protease; N protein: nucleocapsid protein; N7-MTase: guanine-N7 methyl transferase; NCB site: nucleic acid binding site; Nsp: nonstructural protein; PLpro: papain-like protease; RdRp: RNA-dependent RNA polymerase; SAM: S-Adenosylmethionine.

Fig. 8

Molecular models (3D and 2D) of the selected chemical constituents binding to the COVID-19 protein targets. GTA: 7-methyl-GpppA; NCB site: nucleic acid binding site; Nsp16: nonstructural protein16; PLpro: papain-like protease; RdRp with RNA: RNA-dependent RNA polymerase; SAM: S-Adenosylmethionine. Green dotted lines: hydrogen bonds; Red combs: hydrophobic contacts.

According to statistical analysis in the PDB database, hydrophobic contacts and hydrogen bonds were the most frequent type of interactions observed in protein-ligand atomic interactions [29]. As shown in Fig. 8, compound ginsenoside Rf bound to RdRp with RNA formed 8 hydrogen bonds and 11 hydrophobic contacts. Chikusetsusaponin Iva bound to RdRp with RNA formed 6 hydrogen bonds and 7 hydrophobic contacts. Ginsenoside Rh4 bound to PLpro formed 6 hydrogen bonds and 10 hydrophobic contacts. Ginsenoside Rg4 bound to three sites in two molecules: when bound to helicase, it formed 6 hydrogen bonds and 11 hydrophobic contacts; when bound to Nsp16 (GTA site), it formed 8 hydrogen bonds and 10 hydrophobic contacts; when bound to Nsp16 (SAM site), it formed 6 hydrogen bonds and 12 hydrophobic contacts.

The helicase protein is an ideal target to develop anti-viral drugs due to its sequence conservation and indispensability in all coronavirus species [30]. Hence, the helicase-ginsenoside Rg4 complex was selected for further molecular dynamics simulation to verify the stability and dynamic interaction. As shown in Fig. 9 , the root-mean-square deviation curve of the helicase backbone exhibited a smooth fluctuation after 6 ns and maintained an average value of 0.4 nm, which suggested that the combination of helicase and ginsenoside Rg4 is relatively stable. The fluctuation of the amino acid residues represented by the root-mean-square fluctuation was mostly within the range of 0.07 to 0.5 nm. The major fluctuations were observed around residue numbers 346 and 476, and in between 550–600, indicating that the protein had greater residue flexibility in these regions. From the Gibbs free energy landscape, the helicase-ginsenoside Rg4 complex was in its lowest binding free energy (0.0 kJ/mol) in the dark blue region, suggesting a good stability of the ligand–protein complex.

Fig. 9

Molecular dynamics simulation of helicase-ginsenoside Rg4. A: root-mean-square deviation (RMSD) plot; B: root-mean-square fluctuation (RMSF) plot; C: Gibbs free energy landscape. Rg: radius of gyration.

Discussion

Severe COVID-19 patients often develop acute respiratory distress syndrome and multiple organ dysfunction syndromes with a relatively high mortality rate that is related to “cytokine storm” [31]. To date, there are still no specific antiviral drugs that target SARS-COV-2. In the previous multicenter prospective clinical trial, SHG showed a significant therapeutic effect in severe COVID-19 [9]. Moreover, another retrospective study conducted by our group has shown a reduction in the level of inflammatory factors such as IL-6, IL-10 and CRP in the SHG treatment group [32]. It can be inferred that SHG reduces inflammation factors and suppress the inflammatory response in severe COVID-19 patients.

In this study, network pharmacology analysis was used to further investigate the targets and predict the mechanisms of SHG in the treatment of severe COVID-19. A total of 131 potential therapeutic target genes were obtained, including HSP90AA1, MMP9, CXCL8, PTGS2, IFNG, DNMT1, TYMS, MDM2, HDAC3 and ABCB1, which were identified as hub genes. Among them, HSP90AA1 was the most critical target in the PPI network. HSP90AA1 is a highly-conserved molecular chaperone that promotes the maturation, structural maintenance, and regulation of specific target proteins. Recent studies showed that the level of HSP90AA1 was higher in SARS-CoV-2-infected host cells and inhibition of the HSP90 activity reduced both the rate of viral replication and expression of pro-inflammatory cytokines [33]. MMP9 is a secreted protein that plays an essential role in local proteolysis of the extracellular matrix. MMP9 has been widely studied in acute lung injury and reported to have a relationship with COVID-19 [34]. As the upregulation of MMP9 can promote the neutrophils-mediated inflammatory response, MMP9 is considered to be an immune-related gene that can affect the clinical severity of COVID-19 [18]. CXCL8 commonly encodes IL-8, which is a major mediator of the inflammatory response. Cambier et al. [35] reported elevated expression of chemokines, including CXCL1, CXCL8, CXCL11 and CXCL12α, in broncho-alveolar lavage fluids of COVID-19 patients. IFNG encodes a soluble cytokine secreted by both the innate and adaptive immune systems. Elevated levels of IFNG have been found in severe COVID-19 patients and are thought to exacerbate the cytokine storm [36]. PTGS2 encodes the inducible isozyme involved in the prostaglandin biosynthesis associated with mitogenesis and inflammation. Both IFNG and PTGS2 participate in COVID-19 and are associated with the disease severity [37]. We found that most of these hub genes involved in inflammation and immunity have played a crucial role in the development of COVID-19. Hence, these hub genes can be confirmed as therapeutic targets for SHG.

Severe cases of COVID-19 can be defined as sepsis induced by SARS-CoV-2 infection, which causes the maladjusted response of the body and deadly organ failure [38]. In addition, ACE2, a functional receptor for SARS-CoV-2, is widely distributed in multiple organs including the lungs, heart, kidneys, brain, intestine and testes [39]. Although ACE2 is rarely expressed in immune cells, some immune cells can be infected by SARS-CoV-2 because CD147 and CD26, which are strongly expressed in nearly all immune cells, may act as secondary receptors [40], [41]. Therefore, we utilized BioGPS to observe the distribution of all therapeutic targets in the body and constructed the target-organ location network to explore the potential localization of SHG activity in the body. The results revealed that the therapeutic targets are highly related to immune cells, reflecting the immunomodulatory effect of SHG on severe COVID-19. Further, the therapeutic targets are also expressed specifically in the brain, lymphoid tissues, heart, gastrointestinal tract and endocrine tissues, indicating a synergistic therapeutical effect of SHG on the whole body.

GO and KEGG enrichment analyses were carried out to explore the underlying mechanism of SHG in treating severe COVID-19. The BP domain was linked to the responses to various external stimuli, such as drugs, decreased oxygen levels, and regulation of metabolic process. The GO domain suggested that SHG played an important role in the process of host response to and defense against infections. In terms of KEGG analysis, the annotation genes were mainly enriched in neuroactive ligand-receptor interactions and calcium signaling pathways. A previous study has reported that neuroactive ligand-receptor interactions may be involved in the pain response process [42]. Thus, we infer that SHG may have some palliative effect on myalgias, headache and abdominal pain caused by the virus infection. Another study proved that calcium signaling is an integral part of the integrin and chemokine signaling pathways that are vitally involved in synchronizing the neutrophil recruitment process [43]. Pathways like prostate cancer, purine metabolism, alcoholism, gastric cancer, morphine addiction and small cell lung cancer are not closely related to the mechanism of SHG against COVID-19. The immune system-related enrichment pathways, including the IL-17 signaling pathway, TLR signaling pathway, CLR signaling pathway and TCR signaling pathway, are worthy of attention. The IL-17 family is a subset of cytokines that play crucial roles in both acute and chronic inflammatory responses [44]. TLR and CLR function as pattern-recognition receptors that are responsible for detecting microbial pathogens and generating immune responses [45], [46]. TCR-mediated T lymphocyte activation is a key event for effective immune system response. T cell malfunctions in responding to the virus in the early stage of infection are associated with disease progression and poor prognosis [47].

Since the SARS-COV-2 proteins encoded by the viral genome are important for viral replication and transcription [48], the protein binding site is an ideal target to develop anti-viral drugs. To facilitate drug discovery against COVID-19, we screened 10 key chemical constituents with a high degree of connectivity in the herb-constituent-target network and 15 SARS-COV-2 protein targets for molecular docking. The results showed that the key chemical constituents had good affinity towards SARS-COV-2 structural and nonstructural proteins. Compounds including ginsenoside Rf, chikusetsusaponin Iva, ginsenoside Rh4 and ginsenoside Rg4 presented superior affinity interactions with RdRp with RNA, helicase, PLpro and Nsp16, primarily through hydrogen bond and hydrophobic interactions. The molecular dynamics simulation further verified the stability and interaction of helicase and ginsenoside Rg4. These findings demonstrate that SHG may have a broad-spectrum antiviral effect on severe COVID-19 by acting on multiple protein targets. However, the main targets and mechanisms predicted in this research still need to be verified by in vivo and in vitro experiments.

Conclusion

Our study investigated the targets and predicted the mechanisms of SHG action in the treatment of severe COVID-19 using network pharmacology, molecular docking and molecular dynamics simulation. Multiple active chemical constituents in SHG regulated multiple targets to inhibit virus invasion and cytokine storm through multiple signaling pathways, which provides a scientific basis for clinical applications and further experiments. Our research suggests that SHG may be a promising complementary drug for the treatment of severe COVID-19.

Funding

This study was supported by the National Key Research and Development Program (No. 2018YFC1705900), and the Emergency Committee of the World Federation of Chinese Medicine Societies and Shanghai Society of Traditional Chinese Medicine Novel Coronavirus Pneumonia Emergency Tackling Key Project (No. SJZLJZ.N01).

Authors’ contribution

FBJ conceived and designed the research. XXR, ZW, WXX, YHQ, SYT, PYT, WB, PW, SLH and GQ performed the acquisition and analysis of data. XXR analyzed the data. FBJ and ZS supervised the research. XXR and ZW wrote the original draft. FBJ and ZS reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

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.