Network pharmacology and molecular docking analysis on molecular targets and mechanisms of Huashi Baidu formula in the treatment of COVID-19.

PURPOSE: Huashi Baidu formula (HSBDF) was developed to treat the patients with severe COVID-19 in China. The purpose of this study was to explore its active compounds and demonstrate its mechanisms against severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) through network pharmacology and molecular docking.
METHODS: All the components of HSBDF were retrieved from the pharmacology database of TCM system. The genes corresponding to the targets were retrieved using UniProt and GeneCards database. The herb-compound-target network was constructed by Cytoscape. The target protein-protein interaction network was built using STRING database. The core targets of HSBDF were analyzed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The main active compounds of HSBDF were docked with SARS-CoV-2 and angiotensin converting enzyme II (ACE2).
RESULTS: Compound-target network mainly contained 178 compounds and 272 corresponding targets. Key targets contained MAPK3, MAPK8, TP53, CASP3, IL6, TNF, MAPK1, CCL2, PTGS2, etc. There were 522 GO items in GO enrichment analysis (p < .05) and 168 signaling pathways (p < .05) in KEGG, mainly including TNF signaling pathway, PI3K-Akt signaling pathway, NOD-like receptor signaling pathway, MAPK signaling pathway, and HIF-1 signaling pathway. The results of molecular docking showed that baicalein and quercetin were the top two compounds of HSBDF, which had high affinity with ACE2.
CONCLUSION: Baicalein and quercetin in HSBDF may regulate multiple signaling pathways through ACE2, which might play a therapeutic role on COVID-19.

Introduction

Coronavirus disease (COVID-19) was caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a newly discovered coronavirus with a size of 60–140 nm, located in the SARB subgene of the Betacoronavirus family [1]. Moreover, the population's susceptibility to this highly pathogenic coronavirus led to massive global outbreak which turned into an international public health event [2]. There is a late surge of the confirmed cases in European region and region of the Americas. As on 2:00am CEST, 18 April 2020, there have been 2,160,207 confirmed cases of COVID-19, including 146,088 deaths globally [3]. A total of 83,160 cases have been confirmed in China, with 3342 deaths and a case fatality rate of 4% by now. According to an analysis of 44,672 confirmed cases out of 72,314 reported cases by the Chinese Center for Disease Control and Prevention (CDC), the crude case fatality rate of the patients with severe COVID-19 reached 49%. The fatality rate in severe patients with acute respiratory distress syndrome (ARDS) was close to 50%, and even as high as 70% if ARDS reached moderate to severe [4]. The older people or the crowd of underlying medical problems were more likely to infect with severe symptoms. The symptoms of patients were common with fever, shortness of breath, upper airway congestion, dry cough, fatigue, sputum production, myalgia/arthralgia, and breathing difficulties [5]. Critical patients can lead to pneumonia, kidney failure, and even death [5].

No vaccine for SARS-CoV-2 has been published publicly, and there was no medication specific for the treatment of COVID-19 so far [6]. It was proved that Traditional Chinese Medicine (TCM) could obviously shorten fever duration and symptomatic relief of the patients with severe COVID-19 [7-9]. Up to now, ‘three medicines and three formulas’ have been proven effective in treating the COVID-19. Among them, Huashi Baidu formula (HSBDF) was used as an auxiliary medicine for the treatment of patients with severe COVID-19. HSBDF was developed by Chinese Academy of Traditional Chinese Medicine [10]. HSBDF is consisted of glycyrrhiza, apricot kernel, agastache rugosus, magnolia officinalis, atractylodes, amomum tsao-ko, pinellia ternate, poria cocos, rhubarb, astragalus, lepidium seed, red peony root, and ephedra. According to the theory of TCM, the core pathogenesis of COVID-19 was the wet epidemic caused by the cold and humidity outside the lung and spleen, which was transformed into heat and lead to heat stagnation due to the imbalance of qi mechanism and endogenous stagnated heat [11]. Agastache rugosus and ephedra in HSBDF have the functions of detoxifying dampness, clearing heat, and relieving asthma [10]. HSBDF was applied to the treatment of severe patients, and proved to be effective in resisting SARS-CoV-2, eliminating inflammation, and improving immunity [10,12].

However, the mechanism of HSBDF for the treatment of COVID-19 was not clear. HSBDF contained 14 Chinese herbs, and the components of each herb were complex. The therapeutic effect of each herb or each component was not clear.

Network pharmacology was proposed as a promising method to understand herbal formulas [13] and predict potential new drugs or targets for the specific diseases [14-16]. Molecular docking is a significant pathway of structural molecular biology and the computer aided drug design in new medicines [17,18]. Study showed that SARS-CoV-2 binded to angiotensin converting enzyme II (ACE2) receptor with nearly 10–20 times higher affinity than SARS-CoV [17]. The combination of SARS-CoV-2 and ACE2 was the main cause of COVID-19. Therefore, ACE2 and SARS-CoV-2 3CL were regarded as receptors in molecular docking. In this study, we aim to utilize network pharmacology and molecular docking to understand the active compounds of HSBDF, predict their potential targets and signal pathways, and explore the association between the active compounds of HSBDF with ACE2.

Materials and methods

Identification and screening of active compounds

Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, http://tcmspw.com/) was used to authenticate all compounds of the fourteen Chinese medicinal herbs in HSBDF. TCMSP had the ability to identify ingredient–target networks and the large number of herbal items. It was composed of 499 Chinese herbs recorded in the Chinese pharmacopeia with 29,384 compounds and 3311 targets [19]. The names of herbs were used as the key words to retrieve all components. We screened compounds of HSBDF based on absorption, distribution, metabolism, and excretion (ADME), which were strategic processes in drug discovery and development [20]. Oral bioavailability (OB) and drug-likeness (DL) were the most significant pharmacokinetic parameters. In ADME system [21], OB was defined as ‘the degree to which active ingredients are used by the body, including the dose and rate at which they enter the bloodstream’ by Food and Drug Administration (FDA) [22]. The degree of OB largely determined the effect of the compound on disease. DL was used to screen out first-rank compounds and improve candidate compounds in the early period of drug development [14]. In this study, the active compounds in HSBDF were selected according to the criterion of OB ≥30% and DL ≥0.18 [23].

Identification of protein targets

The protein targets associated with active compounds were retrieved from the TCMSP database (http://lsp.nwsuaf.edu.cn/tcmsp.php), which provided information of 6511 drug molecules and 3987 targets as well as the interactions between them [19]. As an authoritative database of protein sequences, UniProt Knowledgebase (UniProtKB) comprised 54,247,468 sequence items [24]. The targets, including the gene names and gene ID, were further extracted using UniProtKB (http://www.uniprot.org).

Construction of component–target gene network

Visual component–target network was established based on aforementioned data sets through Cytoscape 3.7.2 (http://www.cytoscape.org/) to reflect the complex relationships between active compounds and their potential targets [25]. Cytoscape 3.7.2 is an open-source software platform which is used to visualize complicated networks and integrate different types of attribute data [26]. In the network, nodes represent the screened active ingredients and targets, while the connections between the nodes represent the interactions between these biological analyses. The top eight compounds were screened as ligands for molecular docking based on the degree value. The degree value of the molecular represents the number of connections between the molecular and target in the core architecture of the network [27]. The larger the value is, the more likely the component is to become the key ingredient of HSBDF.

Predicting the targets of COVID-19

GeneCards database (https://www.genecards.org/) and Therapeutic Target Database (TTD, https://db.idrblab.org/ttd/) were used to gather the information on COVID-19-associated target genes [28]. GeneCards is a comprehensive database of functions involving proteomics, genomics, and transcriptomics [29]. The keywords ‘novel coronavirus pneumonia,’ ‘cough,’ and ‘fever’ were utilized to screen the COVID-19-associated targets. The names of targets were collected from TTD, which provided information about the therapeutic protein targets, the corresponding ID of each targets and the targeted disease. The common targets of HSBDF and COVID-19 were then gathered as the core targets of HSBDF for COVID-19.

Construction of protein–protein interactions (PPI) network

The core targets of HSBDF were put into STRING (https://string-db.org/cgi/input.pl) to build the PPI network interaction [30]. Cytoscape V3.7.2 [26] was used to construct and visualize the PPI network. CytoNCA, a network topology analysis plug-in in Cytoscape, was used for topological analysis of PPI networks. ‘Degree’ referred to the number of connections of the node in the whole network, which reflected the interaction information between nodes. The value of ‘Degree’ was used as a reference for the importance of the core target.

GO and KEGG pathway enrichment analysis

Webgestalt (www.webgestalt.org) was used to process data and visualize the enrichment results of Gene Ontology (GO) enrichment analysis. GO enrichment analysis included biological processes (BP), molecular functions (MF), and cellular components (CC). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed using the Rpackage ‘ClusterProfifiler’ [31]. Statistical significance threshold of enrichment analysis was set at p ≤ .05.

Component–target molecular docking

The three-dimensional (3D) structure of the SARS-CoV-2 3CL (PDB ID: 6lu7) was downloaded from the RCSB PDB database (https://www.rcsb.org/). The 3D structure of ACE2 (PDB ID: 1r42) was from Qingdao Ocean Science and Technology National Laboratory 2019-nCoV drug target system structure of information sharing platform (http://ncovtarget.qnlm.ac/web/mg/hm). AutoDock Tools 1.5.6 software was used to remove the water molecules, isolate proteins, add the nonpolar hydrogen and calculate Gasteiger charges for the structure and save it as a PDBQT file. PubChem database (https://pubchem.ncbi.nlm.nih.gov/) was used to download the two-dimensional (2D) structures of the top eight compounds. The 2D structure was processed and transformed into PDB format through Chem3D, and they were saved in PDBQT format as docking ligands in AutoDock Tools 1.5.6 software. SARS-CoV-2 3CL protein and ACE2 were used as receptors, and the active compounds were used as ligands. The active site of molecular docking was determined by the ligand coordinate in the target protein complex. The ligand was set to be flexible and the receptor was rigid. Autodock Vina 1.1.2 was used to dock small molecules with ACE2 protein and SARS-CoV-2 3CL protein, respectively. For results of docking, a total of 20 conformations were generated each. The conformation with the best affinity was selected as the final docking conformation and visualized in Pymol 2.3. These active ingredients were compared with lopinavir, ritonavir, and remdesivir.

Results

Active compounds of HSBDG

A total of 222 active compounds were retrieved in TCMSP database. These active compounds were primarily originated from glycyrrhiza (92 compounds), red peony root (29 compounds), ephedra (23 compounds), astragalus (20 compounds), apricot kernel (19 compounds), rhubarb (16 compounds), poria cocos (15 compounds), pinellia ternate (13 compounds), lepidium seed (12 compounds), agastache rugosus (11 compounds), atractylodes (nine compounds), amomum tsao-ko (eight compounds), etc. The main active components of HSBDG are shown in Table 1.

Table 1.

Basic information of some active components of HSBDG.

Mol ID	Chemical component	OB/%	DL	Herb
MOL004990	7,2′,4′-Trihydroxy-5-methoxy-3－arylcoumarin	83.71	0.27	Glycyrrhiza
MOL002311	Glycyrol	90.78	0.67
MOL001925	Paeoniflorin_qt	68.18	0.40	Red peony root
MOL001918	Paeoniflorgenone	87.59	0.37
MOL004328	Naringenin	59.29	0.21	Ephedra
MOL005190	Eriodictyol	71.79	0.24
MOL000378	7-O-methylisomucronulatol	74.69	0.30	Astragalus
MOL000398	Isoflavanone	109.99	0.30
MOL010922	Diisooctyl succinate	31.62	0.23	Apricot kernel
MOL005030	Gondoic acid	30.70	0.20
MOL002293	Sennoside D_qt	61.06	0.61	Rhubarb
MOL000471	Aloe-emodin	83.38	0.24
MOL000282	Ergosta-7,22E-dien-3beta-ol	43.51	0.72	Poria cocos
MOL000300	Dehydroeburicoic acid	44.17	0.83
MOL006967	Beta-D-Ribofuranoside, xanthine-9	44.72	0.21	Pinellia ternate
MOL006957	(3S,6S)-3-(benzyl)-6-(4-hydroxybenzyl)piperazine-2,5-quinone	46.89	0.27
MOL000422	Kaempferol	41.88	0.24	Lepidium seed
MOL000098	Quercetin	46.43	0.28
MOL005911	5-Hydroxy-7,4′-dimethoxyflavanon	51.54	0.27	Agastache rugosus
MOL005890	Pachypodol	75.06	0.40
MOL000186	Stigmasterol 3-O-beta-D-glucopyranoside_qt	43.83	0.76	Atractylodes
MOL000179	2-Hydroxyisoxypropyl-3-hydroxy-7-isopentene-2,3-dihydrobenzofuran-5-carboxylic	45.20	0.20
MOL000096	(−)-Catechin	49.68	0.24	Amomum tsao-ko
MOL000074	(4E,6E)-1,7-bis(4-hydroxyphenyl)hepta-4,6-dien-3-one	67.92	0.24
MOL005980	Neohesperidin	57.44	0.27	Magnolia officinalis
MOL005970	Eucalyptol	60.62	0.32

The construction of herb–compound–target network

After 44 compounds without target were removed, herb–compound–target network contained 463 nodes (including 13 herbs, 178 compounds and 272 genes) and 4081 edges (Figure 1). The larger node meant more importance. According to the degree analysis, the top eight compounds were MOL000098 (quercetin), MOL000422 (kaempferol), MOL000358 (beta-sitosterol), MOL000449 (stigmasterol), MOL000354 (isorhamnetin), MOL002714 (baicalein), MOL004328 (naringenin), and MOL000392 (formononetin), with 886°, 240°, 185°, 120°, 99°, 76°, 72°, and 72°, respectively. More details of these top eight compounds are shown in Table 2.

Figure 1.

Herb-compound-target network of HSBDF (The triangle nodes are composed of all the herbs of HSBDF, which are surrounded with their particular compounds. The octagon nodes represent the compounds of HSBDF. The rhombus nodes, arranged into a rectangular matrix, represent the relative gene targets of HSBDF).

Table 2.

Basic information of the top eight degree of the compounds.

MOL ID	compound	CAS	Molecular formula
MOL000098	Quercetin	117-39-5	2-(3,4-Dihydroxyphenyl)-3,5,7-trihydroxychromen-4-one
MOL000422	Kaempferol	520-18-3	3,5,7-Trihydroxy-2-(4-hydroxyphenyl)chromen-4-one
MOL000358	Beta-sitosterol	474-58-8	Beta-sitosterol 3-O-glucoside_qt
MOL000449	Stigmasterol	83-48-7	(3S,8S,9S,10R,13R,14S,17R)-17-[(2R,5S)-5-ethyl-6-methylhept-3-en-2-yl]-10,13-dimethyl-2,3,4,7,8,9,11,12,14,15,16,17-dodecahydro-1H-cyclopenta[a]phenanthren-3-ol
MOL000354	Isorhamnetin	480-19-3	3,5,7-Trihydroxy-2-(4-hydroxy-3-methoxyphenyl)chromen-4-one
MOL002714	Baicalein	491-67-8	5,6,7-Trihydroxy-2-phenylchromen-4-one
MOL004328	Naringenin	480-41-1	(2S)-5,7-dihydroxy-2-(4-hydroxyphenyl)chroman-4-one
MOL000392	Formononetin	485-72-3	7-Hydroxy-3-(4-methoxyphenyl)-4H-benzopyran-4-one

Prediction results of virus targets and the construction of PPI network

A total of 216 possible targets of COVID-19 were screened through ‘novel coronavirus pneumonia,’ ‘cough,’ and ‘fever.’ These 216 targets were combined with 272 targets of HSBDF to obtain a total of 53 core targets (Figure 2, Table 3). The PPI network showed that there were strong correlations among targets and a complex interlaced network. The network contained a total of 53 nodes, which were the core targets of HSBDF in the treatment of COVID-19. Among them, MAPK3 (45), MAPK8 (44), TNF (44), IL6 (44), and TP53 (44) were considered to be hub genes (Table 3).

Figure 2.

The PPI network of 53 nodes (genes). The larger nodes in the inner ring represent more important hub nodes. The smaller nodes in the outer ring represent the other nodes.

Table 3.

A partial information table for the core targets.

Targets	Degree	Targets	Degree
MAPK3	45.0	MAPK1	42.0
MAPK8	44.0	IL10	41.0
TNF	44.0	CXCL8	39.0
IL6	44.0	PTGS2	38.0
TP53	44.0	CCL2	38.0
CASP3	43.0	IL1B	38.0

GO and KEGG enrichment analyze

To further understand the intersection genes, GO enrichment analysis was conducted. (1) A total 328 top-ranking terms about biological processes were selected, mainly including inflammatory response, response to drug, and cellular response to lipopolysaccharide (Figure 3). The active targets essentially connected to positive regulation of transcription from RNA polymerase II promoter or apoptotic process, including IL6, TNF, and TP53. (2) According to enrichment analysis of cellular components, the targets mainly contained extracellular space, Cytosol, extracellular region, Cytoplasm, Nucleus, and so on. (3) Simultaneously, molecular function terms mainly contained protein binding, identical protein binding, enzyme binding, cytokine activity, protein homodimerization activity, etc.

Figure 3.

The GO enrichment analysis of 53 nodes (BP, MF, CC).

KEGG pathway enrichment was conducted to cluster the major effects that associated to the HSBDF. A total of 20 top-ranking pathway (Figure 4) screened out (p < .05). The main pathway included TNF signaling pathway, PI3K–Akt signaling pathway, pathways in cancer, MAPK signaling pathway as well as some related to cancer.

Figure 4.

Target–pathway network of HSBDF. (The round nodes in the center represent the important hub nodes, and the smaller round nodes in the middle ring represent the other nodes. The outermost triangles represent the related pathways).

Results of molecular docking

Molecular docking was used to verify if the top eight compounds had a significant role in regulation of ACE2. The result showed that all the key compounds in the network had strong affinity with ACE2 protein and SARS-CoV-2 3CL protein. Baicalein (MOL002714) was the most stable active ingredient in SARS-CoV-2 3CL binding. Meanwhile, quercetin (MOL000098) showed strong association with ACE2. The small molecule quercetin showed a compact binding pattern with ACE2 protein active pocket. Quercetin formed four hydrogen bonds with the amino acid residues Lys745, Tyr613, His493, and Asp609, making quercetin and ACE2 form a stable complex (Figure 5). The binding energies of the various compounds are shown in Table 4.

Figure 5.

(a) ACE2 protein-quercetin. (b) SARS-CoV-2 3CL protein-baicalein.

Table 4.

The binding energy of top eight compounds.

CAS	Molecule name	Molecular formula	MW	SARS-CoV-2 3CL docking score (kcal/mol)	ACE2 docking score (kcal/mol)
117-39-5	Quercetin	C15H10O7	302.25	−7.5	−8.4
520-18-3	Kaempferol	C15H10O6	286.25	−7.6	−8.2
83-46-5	Beta-sitosterol	C29H50O	414.79	−7.7	−8.3
83-48-7	Stigmasterol	C29H48O	412.77	−7.0	−8.2
480-19-3	Isorhamnetin	C16H1207	316.28	−7.2	−8.2
491-67-8	Baicalein	C15H10O5	270.25	−7.8	−8.2
480-41-1	Naringenin	C15H12O5	272.27	−7.7	−7.9
485-72-3	Formononetin	C16H12O4	268.28	−7.2	−7.5
192725-17-0	Lopinavir	C37H48N4O5	628.80	−8.0	−8.2
1809249-37-3	Remdesivir	C27H35N6O8P	602.58	−7.6	−7.9
155213-67-5	Ritonavir	C37H48N6O5S2	720.96	−7.8	−8.5

Discussion

In this study, we screened out eight main active components of HSBDF: quercetin, kaempferol, beta-sitosterol, stigmasterol, isorhamnetin, baicalein, naringenin, and formononetin. In molecular docking, these active ingredients were compared with lopinavir, ritonavir, and remdesivir, which were currently recommended for clinical therapeutics [32-34]. Results showed that baicalein showed strong affinity for SARS-CoV-2, and quercetin had strong affinity for ACE2 3CL. It indicated that baicalein and quercetin might play an important role in the treatment of SARS-CoV-2. Baicalein and quercetin were both flavonoids. Flavonoids reduced the barrier dysfunction induced by influenza A virus by inhibiting the NOX4/NF-κB/MLCK pathway, which might be a potential drug for the prevention and treatment of influenza A virus and pulmonary endothelial barrier dysfunction [35]. A study showed that baicalein inhibited the overactivation of the complement system in vivo and improved the acute lung injury induced by influenza a virus [36]. Since both SARS and COVID-19 were caused via binding S-protein to ACE2 [37,38], quercetin acted as a competitive antagonist to inhibit infection of SARS-CoV-2. Other research also reported that quercetin had the functions of reducing capillary brittleness, angiogenesis, detoxification, apoptosis, cell cycle, and antioxidant replication [39].

KEGG enrichment analysis showed that the key targets were mainly concentrated in TNF signaling pathway, PI3K–Akt signaling pathway, MAPK signaling pathway, HIF-1 signaling pathway, NOD-like receptor signaling pathway. The results of PPI network showed that, MAPK3, MAPK8, TNF, IL6, and TP53 were considered to be hub genes. According to these results, we think HSBDF had effect to treat SARS-CoV-2 through the following pathways. (1) TNF signaling pathway: GO enrichment analysis also showed that the major biological processes included inflammatory response in our study. TNF signaling pathway was an important pathway in inflammatory response [40], in which related factor receptors can also induce apoptosis. Quercetin acted as an anti-inflammatory by regulating the TNF signaling pathway [41]. Among the hub genes, IL6 and TNF were involved in TNF signaling pathway. IL6 and TNF were important inflammatory cytokines and commonly involved in the process of inflammation [42]. TNF induced the production of IL-6 and other cytokines, participated in the process of inflammation and oxidative stress [43]. (2) MAPK signaling pathway: Previous studies showed that MAPK signaling pathway participated in the progression of ARDS [44,45]. Many inflammatory factors such as IL-1β, TNF-α, and IL-6 were produced via MAPK signaling pathway [46]. Some anti-inflammatory medications work by targeting MAPK signaling pathway [47,48]. Quercetin was found to regulate the activation of MAPK signaling pathway in retinoblastoma, cardiomyocytes, and chorionic carcinoma cells [49]. Other studies have found that baicalin may inhibit the expression and invasion of cancer cells by inhibiting the p38 MAPK signaling pathway [50]. (3) PI3K–Akt signaling pathway: The PI3K/AKT signaling pathway regulated the activation of inflammatory response cells and the release of inflammatory transmitters to play a role in chronic inflammatory response in the lungs and airways [51]. Quercetin inhibited the PI3K–Akt signaling pathway by inhibiting the expression of AKT1 to silence the anti-apoptotic effect of lung fibroblast, thus realize the treatment of pulmonary fibrosis [52]. Baicalin also played an anti-pulmonary fibrosisrole by the PI3K–Akt pathway [53-55].

There are several limitations in our study. First, our results need to be further verified by experiments. Second, more comprehensive TCM target genes database was needed, which made the results of network pharmacology analysis more reliable. Third, even if the results of network pharmacology and molecular docking were combined, we still could not completely understand the accurate therapeutic mechanism of HSBDF. A comprehensive understanding of HSBDF and COVID-19 depended on the common development of multi-disciplines.

Conclusions

In summary, network pharmacology showed that the main active components of HSBDF, particularly baicalein, and quercetin could act on multiple targets. HSBDF had effect to treat SARS-CoV-2 mainly through the following pathways: TNF signaling pathway, PI3K–Akt signaling pathway, MAPK signaling pathway. Molecular docking showed baicalein and quercetin were the top two compounds which indicated that they might play an important role in the treatment of SARS-CoV-2.