Protection against COVID-19 injury by qingfei paidu decoction via anti-viral, anti-inflammatory activity and metabolic programming.

Qingfei Paidu decoction (QFPD), a multi-component herbal formula, has been widely used to treat COVID-19 in China. However, its active compounds and mechanisms of action are still unknown. Firstly, we divided QFPD into five functional units (FUs) according to the compatibility theory of traditional Chinese medicine. The corresponding common targets of the five FUs were all significantly enriched in Go Ontology (oxidoreductase activity, lipid metabolic process, homeostatic process, etc.), KEGG pathways (steroid biosynthesis, PPAR signaling pathway, adipocytokine signaling pathway, etc.), TTD diseases (chronic inflammatory diseases, asthma, chronic obstructive pulmonary Disease, etc.), miRNA (MIR183), kinase (CDK7) and TF (LXR). QFPD contained 257 specific targets in addition to HCoV, pneumonia and ACE2 co-expression proteins. Then, network topology analysis of the five components-target-pathway-disease networks yielded 67 active ingredients. In addition, ADMET estimations showed that 20 compounds passed the stringent lead-like criteria and in silico drug-likeness test with high gastrointestinal absorption and the median lethal dose (LD50 > 1600 mg/kg). Moreover, 4 specific ingredients (M3, S1, X2 and O2) and 5 common ingredients (MS1, MX16, SX1, WO1 and XO1) of QFPD presented good molecular docking score for 2019-nCov structure and non-structure proteins. Finally, drug perturbation of COVID-19 network robustness showed that all five FUs may protect COVID-19 independently, and target 8 specifically expressed drug-attacked nodes which were related to the bacterial and viral responses, immune system, signaling transduction, etc. In conclusion, our new FUNP analysis showed that QFPD had a protection effect on COVID-19 by regulating a complex molecular network with safety and efficacy. Part of the mechanism was associated with the regulation of anti-viral, anti-inflammatory activity and metabolic programming.

Introduction

2019-novel coronavirus (2019-nCov) outbreak took place in December 2019 and continues to spread around the world. By April 3, 2020, more than 1 million patients have been diagnosed with corona virus disease 2019 (COVID-19) [1]. The virus has a long incubation period, is highly contagious, and is generally susceptible to all types of people, which has a huge negative impact on people's health, economic development, and social stability [2]. However, there is still a lack of effective clinical drugs or vaccine to control the virus.

Traditional Chinese medicine has a good effect on viral infectious pneumonia and has shown a certain effect in the treatment of SARS. On February 7, 2020, the China Health Commission and the Administration of Traditional Chinese Medicine jointly issued a notice recommending formula Qingfei Paidu decoction (QFPD, Herba Ephedrae, Radix Glycyrrhizae, Semen Armeniacae Amarum, Gypsum Fibrosum, Ramulus Cinnamomi, Rhizoma Alismatis, Polyporus Umbellatus, Rhizoma Atractylodis Macrocephalae, Poria, Radix Bupleuri, Radix Scutellariae, Rhizome Pinelliae Preparata, Rhizoma Zingiberis Recens, Radix Asteris, Flos Farfarae, Rhizoma Belamcandae, Herba Asari, Rhizoma Dioscoreae, Fructus Aurantii Immaturus, Pericarpium Citri Reticulatae, Herba Pogostemonis) for the treatment of COVID-19 according to clinical treatment and efficacy. QFPD is a compound prescription in TCM including Ma Xing Shi Gan decoction (MSXG), She Gan Ma Huang decoction (SGMH), Xiao Chai Hu (XCH), and Wu Ling San (WLS), which was first discovered in the classic Treatise on Exogenous Febrile Disease (Shanghan Lun). MXSG (Herba Ephedrae, Radix Glycyrrhizae, Semen Armeniacae Amarum, Gypsum Fibrosum) has been used for the treatment of the common cold, fever, and influenza virus infections via damaging the viral surface structure and inhibiting viral entry [3]. SGMH (Herba Ephedrae, Rhizome Pinelliae Preparata, Rhizoma Zingiberis Recens, Radix Asteris, Flos Farfarae, Rhizoma Belamcandae, Herba Asari) is a classical prescription for the treatment of flu-like symptoms, asthma, inflammation, tonsillitis and sore throat [4]. XCH (Radix Glycyrrhizae, Radix Bupleuri, Radix Scutellariae, Rhizome Pinelliae Preparata, Rhizoma Zingiberis Recens) possesses antiviral [5] and various anticarcinogenic properties [6]. WLS (Ramulus Cinnamomi, Rhizoma Alismatis, Polyporus Umbellatus, Rhizoma Atractylodis Macrocephalae, Poria), a famous Chinese prescription for nephritic syndrome, can improve kidney excretion function and inhibit inflammatory response [7]. These researches indicate that MXSG, SGMH, XCH and WLS may be functional units of formula QFPD. Previous studies have focused on the mechanism of compound prescription based on a single traditional Chinese medicine. However, it may not reflect functional compatibility mechanism of traditional Chinese medicine. Therefore, it is worthy of comparing the similarities and differences of different QFPD functional units in the treatment of COVID-19, including MXSG, SGMH, XCH, WLS and Others.

QFPD contains a total of 21 traditional Chinese medicines, and it is difficult to elucidate the complex mechanism of QFPD on COVID-19 by traditional pharmacological methods due to the multi-components and multi-targets of the formula. Network pharmacology, a new method in recent years, can integrate interactions of drugs, targets, pathways and diseases into a biological network system [8]. Therefore, more and more TCM researchers have begun to use network pharmacology to explore the material basis of TCM, and to reveal the overall comprehensive effects of multi-path, multi-component and multi-target of TCM prescription and its treatment of diseases [9,10]. More importantly, previous study reported that disease conditions can be more fragile than health systems against various perturbations for the un-optimized system [11]. So the formula may be more effective for COVID-19 disease via the stronger effects on the reduction of the robustness of the COVID-19 disease network [12].

In our study, since MSXG, SGMH, XCH and WLS have been independently used for the treatment of viral infectious pneumonia, this study firstly screened out major effective compounds from five functional units respectively. Then we offered a new understanding of the functional units mechanism of QFPD against COVID-19 by a novel functional units of network pharmacology (FUNP) approach and formula perturbation analysis, and provided a combination strategy to explore mechanisms of inter-ingredients interactions from a holistic perspective.

Materials and methods

Data preparation

Compounds of the main herb in formula MSXG, SGMH, XCH, WLS and Others were searched in TCMSP [13], and screened based on drug-likeness (DL) ≥0.18 [14] and oral bioavailability (OB) ≥30 % [15]. Then, the corresponding Pubchem CIDs of the compounds were retrieved from the Pubchem database [16]. Finally, BATMAN-TCM [17], an bioinformatics analysis tool for studying TCM’s molecular mechanisms, was used to identify potential target genes of the active components (uploaded by Pubchem CIDs). To make the results more credible, we set the cutoff score ≥ 30 as the standard. Finally, to discovery the co-differentially presented targets in the five formulae, we conducted pan-formula analysis using Venn diagrams (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Functional and pathway enrichment analyses of QFPD targets

To better understand the functional involvements of MSXG, SGMH, XCH, WLS and Others targets, bioinformatics analyses of multiple formulae targets were first performed, including Gene Ontology (GO) function term, KEGG biological pathway and OMIM/TTD disease enrichment analyses. Then, kinase, microRNA and transcriptional factor (TF) enrichment analyses of the five formulae targets were conducted using the tool WebGestalt (http://bioinfo.vanderbilt.edu/webgestalt) [18] and the bubble and chord plot map were drawn with the R language ggplot2 and GOplot installation package. P-values were adjusted for multiple testing by Benjamini-Hochberg adjustment.

Construction of PPI network and MCODE modules analysis

To further explore the pharmacological mechanisms, five PPI networks were built including: MSXG, SGMH, XCH, WLS and Others targets PPI network. Specifically, the five kinds of target proteins were respectively uploaded to Metascape to build PPI networks, with the species limited to “Homo sapiens”. Next, MCODE analysis [19], a method for finding densely connected modules in PPI networks, was carried out by Cytoscape 3.2.1 (http://www.cytoscape.org/) [20]. Finally, KEGG (Kyoto Encyclopedia of Genes and Genomes) signaling pathway enrichment analysis was further conducted on the identified functional modules of MSXG, SGMH, XCH, WLS and Others targets PPI networks, respectively.

Network construction

Based on the five formulae’s active components, BATMAN-TCM was used to set up five networks of components-target-pathway-disease (MSXG, SGMH, XCH, WLS and Others). To emphasize the important elements of the five networks, we only exhibited the hub targets according to the default criteria (targets with no fewer than 6, 5, 8, 7 and 4 linking compounds for MSXG, SGMH, XCH, WLS and Others, respectively). Finally, these important linking compounds of MSXG, SGMH, XCH, WLS and Others networks were obtained for further analysis.

ADMET evaluation of the predicted active compounds

Based on the SwissADME database [21], the physicochemical properties of the active components was predicted, including molecular weight (MW), rotatable bonds count, H-bond acceptors and donors count, TPSA and leadlikeness violations. Second, pharmacokinetic properties was predicted through pkCSM database [22], which contained the absorption (Caco-2 cell permeability, HIA and skin permeability), distribution (VDss, unbound fraction, blood-brain barrier and central nervous system permeability), excretion (total clearance and renal OCT2 substrate) and toxicity (AMES toxicity, maximum tolerated dose, hERG I inhibitor, hERG II inhibitor, oral rat acute toxicity (LD50), hepatotoxicity, skin sensitisation, and minnow toxicity).

Molecular docking

To facilitate drug discovery against COVID-19, we used COVID-19 Docking Server (https://ncov.schanglab.org.cn/index.php) [23] to predict the binding modes between 12 COVID-19 targets and the 20 lead-likeness of QFPD. Specifically, the 10 nonstructural and 2 structural proteins of 2019-nCov were collected (Mpro, PLpro, nsp12 [RdRp with RNA], nsp12 [RdRp without RNA], nsp13 [Helicase ADP site], nsp13 [Helicase NCB site], nsp14 [ExoN], nsp14 [N7-MTase], nsp15 [endoribonuclease], nsp16 [2′-O-MTase], N protein NCB site and E protein [ion channel]); and the corresponding Protein Data Bank (PDB)codes were 6LU7, 4OW0, 3H5Y (with RNA), 3H5Y (without RNA), 6JYT (ADP site), 6JYT (NCB site), 5C8S (ExoN),5C8S (N7-MTase), 2RHB, 2XYR, 4KYJ, and 5 × 29, respectively. Finally, Discovery Studio software elucidated the 14 best docking results between compounds and the COVID-19 target proteins.

ACE2 and CD147 expression across tissues and co-expression genes

To understand the expression and distribution of ACE2 and CD147 across tissues, a radar plot including 53 tissues was performed through COXPRESdb [24]. And the top 200 co-expression genes of ACE2 and CD147 (P < 1E-16) were obtained, respectively. Then, text mining method from the literature was used to screen for pneumonia-associated genes through COREMINE (http://www.coremine.com/). In addition, co-expression genes of ACE2 in colonic epithelial cells [25] and HCoV-associated host proteins with references [26] were obtained. Finally, we performed UpsetView analysis (http://www.ehbio.com/ImageGP/) between these five sets of proteins and QFPD targets.

Validation of drug positioning for QFPD against COVID-19 via the robustness of disease network

Since QFPD effects on COVID-19 via multi-component and multi-target, we evaluate the potential efficacy of QFPD through TCMATCOV platform, which uses the quantitative evaluation algorithm of multi-target drugs to disturb the disease network. Specifically, the disturbing effect of drugs on diseases is simulated by deleting disease network nodes. The disturbance rate of drugs is calculated by comparing the changes of network topology characteristics before and after drug intervention, which is used to evaluate the intervention effect of drugs on diseases. Firstly, COVID-19 disease network was constructed based on specific cytokines of COVID-19 [27] and differentially expressed genes of SARS (GSE36969, GSE51387, GSE68820). Then, this platform uses four kinds of network topology characteristics to evaluate the robustness of COVID-19 network, including network average connectivity, network average shortest path, connectivity centrality and compactness centrality. And the five formulae (MSXG, SGMH, XCH, WLS and Others) disturbance scores are calculated according to the changes before and after drug intervention. Finally, the disturbance effect of the five formulae on the COVID-19 network was compared with null models with the total score of the disturbance, and the higher the value is, the higher the damage degree of drugs to the stability of the network is [12]. We take Banxia tianma baizhu decoction (BXTM) as negative control; and another efficient formula Yi du bi fei decoction (YDBF) as positive control.

Result

Prediction of active components and potential targets of QFPD

Firstly, the DL ≥ 0.18 and OB ≥ 30 %.s were set as the standard to screen the chemical components obtained through online database TCMSP. Specifically, a total of 175 effective components of QFPD were screened from the TCM database, including 82 species of MSXG, 35 species of SGMH, 105 species of XCH, 21 species of WLS and 32 species of Others (Table 1 ). Among these effective components, 89 (50.86 %) components exited in more than two formulae; CID5280343 and CID5280794 were owned by MXSG, Others, SGMH and XCH; CID12303645 was owned by MXSG, Others, WLS and XCH (Fig. 1 A). Secondly, a total of 300 targets of QFPD were screened from the BATMAN-TCM database, including 192 targets of MSXG, 201 targets of SGMH, 221 targets of XCH, 96 targets of WLS and 99 targets of Others. Among these proteins, 21 (7%) targets exited in five formulae (Fig. 1B).

Table 1

Effective components of QFPD.

Formula	N	PubChem_Cid
MSXG	82	CID10090416,CID10542808,CID10881804,CID11267805,CID114829,CID11558452,CID11602329,CID11975273,CID120074,CID12303645,CID124049,CID124052,CID13965473,CID14604077,CID14604078,CID14604081,CID15228663,CID15380912,CID162412,CID177149,CID193679,CID197678,CID23724664,CID25015742,CID268208,CID336327,CID354368,CID3764,CID439246,CID440833,CID442411,CID44257530,CID480774,CID480780,CID480787,CID480859,CID480873,CID49856081,CID503731,CID503737,CID5280343,CID5280378,CID5280448,CID5280544,CID5280794,CID5280863,CID5281619,CID5281654,CID5281789,CID5282768,CID5282805,CID5312521,CID5316900,CID5317300,CID5317478,CID5317479,CID5317480,CID5317652,CID5317768,CID5317777,CID5318437,CID5318585,CID5318679,CID5318869,CID5318998,CID5318999,CID5319013,CID5320083,CID5460988,CID5481234,CID5481948,CID5481949,CID5997,CID636883,CID637112,CID64971,CID6918970,CID73205,CID928837,CID9927807,CID15840593,CID15228662
SGMH	35	CID1135,CID1174,CID3026,CID5789,CID6782,CID6998,CID8437,CID8679,CID13625,CID117158,CID159225,CID185,CID10019512,CID11438306,CID11869417,CID11870462,CID12315507,CID162350,CID16726037,CID222284,CID389888,CID3902,CID440833,CID5280343,CID5280445,CID5280544,CID5280794,CID5280863,CID5281331,CID5281605,CID5281616,CID5281628,CID5281654,CID5281779,CID5282768,CID5315890,CID5316876,CID5320945,CID5484202,CID5488781,CID5491637,CID64959,CID64982,CID676152,CID71307581,CID72307,CID13688752
XCH	105	CID10090416,CID10542808,CID10881804,CID11267805,CID11438306,CID114829,CID11558452,CID11602329,CID117443,CID12303645,CID124049,CID124052,CID124211,CID13965473,CID14135323,CID14604077,CID14604078,CID14604081,CID15228662,CID15228663,CID15380912,CID156992,CID158311,CID15840593,CID159029,CID161271,CID162412,CID177149,CID185034,CID193679,CID197678,CID222284,CID23724664,CID25015742,CID25721350,CID268208,CID336327,CID354368,CID373261,CID3764,CID389001,CID389888,CID439246,CID442411,CID44257530,CID44258628,CID480774,CID480780,CID480787,CID480859,CID480873,CID49856081,CID503731,CID503737,CID5280343,CID5280378,CID5280442,CID5280448,CID5280794,CID5280863,CID5281605,CID5281619,CID5281654,CID5281674,CID5281703,CID5281789,CID5282768,CID5312521,CID5316900,CID5317300,CID5317478,CID5317479,CID5317480,CID5317652,CID5317768,CID5317777,CID5318437,CID5318585,CID5318679,CID5318869,CID5318998,CID5318999,CID5319013,CID5319042,CID5319252,CID5320083,CID5320315,CID5320399,CID5321865,CID5322059,CID5322078,CID5460988,CID5481234,CID5481948,CID5481949,CID5484202,CID636883,CID637112,CID64959,CID64971,CID64982,CID73205,CID821279,CID928837,CID9927807
WLS	21	CID10181133,CID10743008,CID12303645,CID14036811,CID15225964,CID15226717,CID15976101,CID182232,CID222284,CID44575602,CID5283628,CID5471851,CID5471852,CID56668247,CID6436630,CID712316,CID73402,CID9064,CID9805290,CID14448075,CID14236575
Others	32	CID10212,CID11824478,CID122159,CID12303645,CID14057197,CID145659,CID17897,CID33934,CID373261,CID40429858,CID42607889,CID439246,CID442834,CID443024,CID5280343,CID5280445,CID5280794,CID5281326,CID5281617,CID5281781,CID5319406,CID5320621,CID5495928,CID5997,CID631170,CID632135,CID676152,CID712316,CID72344,CID79730,CID1149877,CID45359875

Fig. 1

Venn diagram of the five formulae’ active compounds and targets. A: compounds, B: targets.

As shown in Fig. 2 , the 11 enriched GO terms of the targets in all five formulae were found, such as oxidoreductase activity, lipid metabolic process, lipid binding, small molecule metabolic process, homeostatic process, signal transducer activity and cell proliferation. Furthermore, the results of pathway enrichment analysis showed that the 7 KEGG pathways were significantly related to more than 4 formula groups, including dteroid biosynthesis, adipocytokine signaling pathway, neuroactive ligand-receptor interaction, steroid hormone biosynthesis, PPAR signaling pathway, arginine and proline metabolism and ABC transporters (Fig. 3 A). In addition, TTD analysis showed that the 8 diseases were significantly association with more than 3 formula groups, such as acne, Behcet'S disease, benign prostate hyperplasia, intrahepatic cholestasis and brain injury (Fig. 3B). However, the five formulae contained their specific (MSXG, SGMH, XCH, WLS and Others) GO, KEGG and TTD terms. For example, neurological system process and beta-Alanine metabolism terms were specific for MXSG; membrane organization and parasitic infections of the eye terms for SGMH; circadian rhythm and hypertension terms for XCH; nucleic acid binding transcription factor activity, ovarian steroidogenesis and chronic inflammatory diseases terms for WLS; fat digestion and absorption terms for Others.

Fig. 2

Bubble plot of the GO analysis of the five formulae’ targets.

Fig. 3

Bubble plot of the KEGG/TTD analysis of the five formulae’ targets. A: KEGG, B: TTD.

In the prediction of miRNAs in QFPD targets, MIR-183 and MIR-130A/B/301 were the highest linking terms to bind the five formulae targets and formulae Others was the highest group to bind miRNAs (Fig. 4 A). In addition, kinase prediction revealed CDK7 were significantly enriched in formulae MSXG, SGMH, XCH, WLS and Others (Fig. 4B). Finally, TF analysis showed that LXR was the highest linking TF to bind the four formulae targets and formulae WLS was the highest group to bind TFs (Fig. 4C).

Fig. 4

The miRNA, kinase and TF analysis of the five formulae’ targets by WebGestalt. Chord plot showing the five formulae’ targets present in the represented enriched miRNA, kinase and TF terms. Outer ring shows miRNA/kinase/TF term and log2 enrichment ratio (left) or five formulae grouping (right). Chords connect miRNA/kinase/TF term with formulae groups. A: miRNA, B: kinase, C: TF.

To further explore the functional relationship among five formulae, PPI networks were constructed through Metascape, and visual composition carried out by Cytoscape. Firstly, the potential 192 target genes of MXSG were analyzed by PPI network, and the results showed that there were 144 nodes and 510 edges, which represented the interaction between protein and function. The MXSG PPI network function module was confirmed by the MCODE plug-in and a list of the corresponding meaningful modules presented (Fig. 5 A). 3 modules scores were > 2.5. Module 1 (score: 5.769) consisted of 13 nodes and the seed gene was COX7A1; Module 2 (score: 4.429) consisted of 14 nodes and the seed gene was ALDH1A1; module 3 (score: 5.0) consisted of 11 nodes and the seed gene was CNR2. KEGG pathway enrichment analysis showed that MXSG modules were enriched in neuroactive ligand-receptor interaction, calcium signaling pathway, inflammatory mediator regulation of TRP channels, et.al.

Fig. 5

KEGG analysis of MCODE modules. MCODE analysis was performed after the construction of the five formulae’ targets PPI; then, KEGG analysis was conducted on the MCODE modules. A: MXSG, B: Others, C: WLS, D: SGMH, E: XCH.

Secondly, the potential 99 target genes of Others were analyzed by PPI network, and the results showed that there were 77 nodes and 194 edges. Only 1 module score were > 2.5 (Fig. 5B). Module 1 (score: 5.769) consisted of 13 nodes and the seed gene was COX7A1. KEGG pathway enrichment analysis showed that Others modules were enriched in huntington's disease, glycolysis / gluconeogenesis, Notch signaling pathway, et.al.

Thirdly, the potential 96 target genes of WLS were analyzed by PPI network, and the results showed that there were 60 nodes and 143 edges. Only 1 modules score were > 2.5 (Fig. 5C). Module 1 (score: 2.706) consisted of 17 nodes and the seed gene was CNR2. KEGG pathway enrichment analysis showed that WLS modules were enriched in thyroid hormone signaling pathway, adipocytokine signaling pathway, neuroactive ligand-receptor interaction, et.al.

Fourthly, the potential 201 target genes of SGMH were analyzed by PPI network, and the results showed that there were 153 nodes and 505 edges. 3 modules scores were > 2.5 (Fig. 5D). Module 1 (score: 3.529) consisted of 17 nodes and the seed gene was ACSS1; Module 2 (score: 4.5) consisted of 7 nodes and the seed gene was CNR2; module 3 (score: 3.5) consisted of 8 nodes and the seed gene was PRKCG. KEGG pathway enrichment analysis showed that SGMH modules were enriched in insulin resistance, adipocytokine signaling pathway, Th17 cell differentiation, et.al.

At last, the potential 221 target genes of XCH were analyzed by PPI network, and the results showed that there were 166 nodes and 643 edges. 5 modules scores were > 2.5 (Fig. 5E). Module 1 (score: 5.769) consisted of 13 nodes and the seed gene was COX7A1; Module 2 (score: 2.769) consisted of 13 nodes and the seed gene was RRM1; module 3 (score: 5.5) consisted of 12 nodes and the seed gene was CNR2; module 4 (score: 2.909) consisted of 11 nodes and the seed gene was ACSS1; module 5 (score: 3.0) consisted of 7 nodes and the seed gene was FFAR1. KEGG pathway enrichment analysis showed that XCH modules were enriched in calcium signaling pathway, cGMP-PKG signaling pathway, neuroactive ligand-receptor interaction, et.al.

After using the BATMAN-TCM, we constructed five ingredients-target-pathway-disease networks, including MSXG, SGMH, XCH, WLS and Others. In order to emphasize the important network elements, we showed the networks that exhibit those targets with larger than 6, 5, 8, 7 and 4 linking compounds for MSXG, SGMH, XCH, WLS and Others, respectively (Fig. 6 ). MSXG network contained 31 key components, 50 proteins and 17 pathways; SGMH network contained 15 key components, 20 proteins and 8 pathways; WLS network contained 18 key components, 9 proteins and 4 pathways; XCH network contained 32 key components, 15 proteins and 12 pathways; Others network contained 10 key components, 13 proteins and 3 pathways. To find the potential drugs of formulae QFPD for COVID-19, a total of 67 hub components were used for ADMET analysis.

Fig. 6

The component-target-pathway-disease network. Purple polygons: PubChem ID of QFPD compounds; blue pentagrams: QFPD targets; yellow circles: KEGG pathway; red square: Therapeutic Target Database (TTD) disease term, green square: Online Mendelian Inheritance in Man (OMIN) disease term. A: MXSG, B: SGMH, C: WLS, D: XCH, E: Others. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

ADMET evaluation of the 67 key compounds

Since in silico ADMET prediction can help early drug design and evaluation, ADMET properties of the 67 key compounds were predicted by SwissADME and pkCSM. Chemical properties including molecular weight (MW), rotatable bonds count, H-bond acceptors and donors count, TPSA and leadlikeness violations were calculated by SwissADME and shown as Fig. 8A. It is worth mentioning that 21 (31.34 %) compounds passed the stringent lead-like criteria (250 g/mol ≤ MW ≤ 350 g/mol, XLOGP ≤ 3.5 and rotatable bonds ≤ 7), which are excellent candidates for drug discovery (Fig. 7 A). And these lead-likeness compounds were further predicted by pkCSM, with the exception of S3 (low gastrointestinal absorption)

Fig. 8

Schematic (3D and 2D) representation that molecular model of specific compounds of each formulae with COVID-19 proteins. A: M3 and E protein [ion channel], B: M3 and nsp13 [Helicase NCB site], C: S1 and nsp13 [Helicase ADP site], D: S1 and PLpro, E: X2 and Mpro, F: O2 and Mpro. M: MXSG, S: SGMH, X: XCH, O: Others.

Fig. 7

Chemical properties statistics of hub components in the formulae. A: Molecular weight, B: rotatable bond count, C: H-bond acceptors count, D: H-bond donors count, E: topological polar surface area (TPSA), F: leadlikeness violations, G: pharmacokinetic and toxicity evaluated parameters of 20 leadlikeness compounds by pkCSM; green = good, yellow = tolerable, red = bad. Caco2: Caco-2 Permeability,HIA: Intestinal Absorption (Human), Skin: Skin Permeability, VDss: volume of distribution, FU: Fraction Unbound (Human), BBB: Blood Brain Barrier permeability, CNS: Central Nervous System permeability,TC: Total Clearance, OCT2: Renal Organic Cation Transporter 2, AMES: AMES toxicity, MTDD: Maximum Tolerated Dose (Human), hERG I/II: hERG I and II Inhibitors, LD50: Oral Rat Acute Toxicity (LD50), HT: Hepatotoxicity, SS: Skin Sensitisation, MT: Minnow toxicity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Regarding the absorption parameters, all 20 compounds (Table 2 ) presented a promising oral availability including the optimal Caco-2 cell permeability, HIA and skin permeability. The drug distribution results showed that most of the compounds distributed in tissue (VDss> 0.45: tissue, VDss <−0.15: plasma) with good unbound fraction scores, thus becoming available to interact with the pharmacological target. Only compound W5 and W11 were entirely unable to penetrate the blood-brain barrier (BBB) and central nervous system (CNS). In addition, 15 compounds presented a good renal elimination and were not substrates of the renal organic cation transporter 2 (OCT2). Finally, 14 compounds did not present any particular toxicity problems including AMES toxicity, maximum tolerated dose, hERG I inhibitor, hERG II inhibitor, oral rat acute toxicity (LD50), hepatotoxicity, skin sensitisation, and minnow toxicity (Fig. 7B).

Table 2

20 potential active compounds from QFPD.

Pubchem	Molecular	Name	Structure	Pubchem	Molecular	Name	Structure
CID6918970	M3	ZINC5356864		CID10019512	S5	3-O-Methylviolanone
CID336327	M5	Medicarpin		CID9064	W5	Cianidanol
CID14057197	O1	–		CID182232	W11	(+)-Epicatechin
CID42607889	O2	Alysifolinone		CID25721350	X1	ZINC13130930
CID3902	S1	letrozole		CID14135323	X2	(2S)-dihydrobaicalein
CID821279	X4	ZINC338038		CID439246	MXO1	naringenin
CID440833	MS1	Leucocyanidol		CID676152	SO1	SR-01,000,767,148
CID177149	MX16	(+)-Vestitol		CID11438306	SX1	cyclo(L-Tyr-l-Phe)
CID114829	MX17	Liquiritigenin		CID712316	WO1	(-)-taxifolin
CID928837	MX8	ZINC519174		CID373261	XO1	Eriodyctiol (flavanone)

M: MXSG, S: SGMH, X: XCH, O: Others.

The application of COVID-19 docking server and Discovery Studio software elucidated the interactions between the 20 lead-likeness compounds (S1, W5, MX17, MX16, W11, M5, XO1, MXO1, SO1, WO1, X4, MX8, M3, S5, SX1, O1, X2, X1, O2, MS1) and the 10 nonstructural and 2 structural proteins (Mpro, PLpro, nsp12 [RdRp with RNA], nsp12 [RdRp without RNA], nsp13 [Helicase ADP site], nsp13 [Helicase NCB site], nsp14 [ExoN], nsp14 [N7-MTase], nsp15 [endoribonuclease], nsp16 [2′-O-MTase], N protein NCB site, E protein [ion channel]). The docking scores were depicted in Table 3, Table 4 . The smaller of docking score, the lower of energy would be required, which means the binding between the compounds and the targets are stronger. There are 9 compounds presenting better bonding ability than other compounds.

Table 3

Docking score between specific ingredients of QFPD and 2019-nCov proteins.

Molecule	M3	M5	S1	S5	W11	W5	X1	X2	X4	O1	O2
Main Protease	−7.7	−7.3	−6.8	−6.9	−7.2	−7.5	−7.1	−7.9	−7.3	−7.1	−7.4
Papain-like protease	−8.7	−8.5	−9.9	−7.7	−8	−8	−8.4	−8.7	−8.2	−8.6	−8.2
RdRp with RNA	−8.4	−8.5	−8.5	−8.3	−9.1	−9.1	−8.4	−8.2	−8.2	−8.3	−8.6
RdRp without RNA	−6.8	−6.8	−6.9	−6.8	−6.7	−6.9	−7.1	−7.2	−6.5	−7.3	−7.1
Helicase ADP site	−6.3	−6.3	−7.3	−6.5	−6	−6	−6.3	−6.5	−6.2	−6.1	−6.5
Helicase NCB site	−7.9	−6.9	−7.2	−7.1	−7.2	−7.2	−7.4	−7.5	−7.3	−7.4	−7.4
Nsp14(ExoN)	−6.8	−6.6	−6.3	−6.4	−6.9	−6.8	−6.7	−6.9	−6.4	−6.6	−6.9
Nsp14(N7-MTase)	−8.8	−8.1	−8.5	−7.5	−8.4	−8.4	−8.3	−8.5	−8	−8.3	−8.5
Nsp15(endoribonuclease)	−6.6	−6.3	−6.3	−5.9	−6.2	−6.2	−6.2	−6.3	−6.2	−6.4	−6.4
Nsp16(2′-O-MTase)	−7.5	−7.4	−7.5	−7.2	−8.2	−8.2	−7.7	−7.9	−7.7	−7.9	−8.4
N protein NCB site	−7.6	−7.5	−7.8	−7.6	−7.6	−7.6	−8	−8	−7.5	−7.6	−7.6
E protein(ion channel)	−8.1	−7	−7.8	−6.7	−6.4	−6.4	−7.2	−7.3	−7.2	−7.1	−6.8

M: MXSG, S: SGMH, X: XCH, O: Others.

Table 4

Evaluation of the effect of QFPD on the robustness disturbance of COVID-19 network.

Topology	MXSG	SGMH	XCH	WLS	Others	BXTM	YDBF
DTS	25.66	26.71	21.02	17.64	23.16	14.52	22.71
AC	−4.63	−5.21	−3.02	−3.49	−5.38	−2.32	−3.78
APL	13.35	13.40	9.96	6.49	10.66	4.59	11.17
CoC	−1.44	−1.64	−1.59	−1.30	−1.23	−1.15	−1.25
ClC	−6.24	−6.46	−6.45	−6.36	−5.88	−6.46	−6.52

Average connectivity (AC), Connection centrality (CoC), Closeness centrality (ClC): the larger the quotient is, the more stable the network and, the less the influence made by the drug. Disturbance total score (DTS), Average length of shortest path (APL), : the larger the quotient is, the less stable the network and, the larger the influence made by the drug. Negative control formula: BXTM. Positive control formula: YDBF.

M3 (Fig. 8A), a specific compound in formulae MXSG, showed eight interactions with E protein [ion channel] including Pi-sigma, Pi-alkyl and Alkyl, which were connected with TYR 57, ALA 32, ILE 46 and PRO 54, etc.; additionally, M3 (Fig. 8B) showed five interactions with nsp13 [Helicase NCB site] including Unfavorable Donor-Donor, Pi-alkyl and Alkyl, which were connected with ASN 559, ARG 409, LEU 42 and PRO 406. S1 (Fig. 8C), a specific compound in formulae SGMH, showed seven interactions with nsp13 [Helicase ADP site] including H-bond interactions, van der waals, Amide-Pi stacked and Pi-alkyl, which were connected with ALA 313, ASP 374, GLN 537 and SER 289, etc.; additionally, S1 (Fig. 8D) showed five interactions with PLpro including Pi-anion, Pi-Pi stacked, Pi-Pi T-shaped and Pi-alkyl, which were connected with TYR 264, ASP 164, TYR 268 and PRO 248. X2 (Fig. 8E), a specific compound in formulae XCH, showed seven interactions with Mpro including H-bond interactions, Pi-Donor hydrogen bond and Pi-alkyl, which were connected with MET 165, GLU 166, LEU 141 and CYS 145, etc. O2 (Fig. 8F), a specific compound in formulae Others, showed seven interactions with Mpro including H-bond interactions, Carbon hydrogen bond, Pi-anion, Pi-sulfur and Pi-alkyl, which were connected with MET 131, GLY 71, LEU 100 and CYS 115, etc.

MS1 (Fig. 9 A), a compound in formulae MXSG and SGMH, showed eleven interactions with N protein NCB site including H-bond interactions, Pi-Donor hydrogen bond, Pi-sigma, Pi-Pi stacked and Pi-alkyl, which were connected with SER 51, THR 109, ALA 50 and PRO 42, etc.; additionally, MS1 (Fig. 9B) showed five interactions with nsp14 [ExoN] including H-bond interactions and Pi-Pi stacked, which were connected with GLU 92, PHE 190, ASP 273 and VAL 91, etc. MX16 (Fig. 9C), a compound in formulae MXSG and XCH, showed seven interactions with nsp15 [endoribonuclease] including H-bond interactions, Alkyl and Pi-alkyl, which were connected with PRO 343, VAL 275, LYS 344 and SER 293, etc. SX1 (Fig. 9D), a compound in formulae SGMH and XCH, showed two interactions with nsp14 [N7-MTase] including Pi-Pi stacked and Pi-alkyl, which were connected with PHE 426; additionally, SX1 (Fig. 9E) showed five interactions with nsp15 [endoribonuclease] including H-bond interactions, Alkyl and Pi-alkyl, which were connected with LYS 344, LYS 289, VAL 291 and PRO 343. WO1 (Fig. 9F), a compound in formulae WLS and Others, showed seven interactions with nsp16 [2′-O-MTase] including H-bond interactions, Carbon hydrogen bond, Pi-Pi T-shaped, Pi-alkyl and Pi-anion, which were connected with PHE 149, CYS 115, ASP 99 and SER 74, etc.; additionally, WO1 (Fig. 9G) showed seven interactions with nsp12 [RdRp without RNA] including H-bond interactions, Carbon hydrogen bond, Unfavorable Donor-Donor, Pi-cation and Pi-anion, which were connected with THR 556, ARG 553, ASP 623 and SER 682, etc. XO1 (Fig. 9H), a compound in formulae XCH and Others, showed ten interactions with nsp12 [RdRp with RNA] including H-bond interactions, Pi-Donor hydrogen bond, Pi-Pi T-shaped and Pi-alkyl, which were connected with CYS 813, GLY 590, LYS 593 and ASP 865, etc.

Fig. 9

Schematic (3D and 2D) representation that molecular model of common compounds of the five formulae with COVID-19 proteins. A: MS1 and N protein NCB site, B: MS1 and nsp14 [ExoN], C: MX16 and nsp15 [endoribonuclease], D: SX1 and nsp14 [N7-MTase], E: SX1 and nsp15 [endoribonuclease], F: WO1 and nsp16 [2′-O-MTase], G: WO1 and nsp12 [RdRp without RNA], H: XO1 and nsp12 [RdRp with RNA]. MS: MXSG and SGMH, MX: MXSG and XCH, SX: SGMH and XCH, WO: WLS and Others, XO: XCH and Others.

Since 2019-nCov may enter other tissues and organs through ACE2 and CD147 binding, we firstly explored the expression and distribution of ACE2 and CD147 across 53 tissues. Fig. 10 A showed that the 5 top expression tissues of ACE2 were terminal ileum, testis, visceral (omentum), left ventricle and kidney cortex, which are 3 fold change higher than lung. And the 5 top expression tissues of CD147 were testis, visceral (omentum), left ventricle, aorta, atrial appendage and transformed fibroblasts. Then, to further understand whether QFPD only targets pneumonia or 2019-nCov, we obtained 200 co-expression genes of ACE2, 200 co-expression genes of CD147, 470 pneumonia-associated proteins, 119 HCoV-associated host proteins, and 476 co-expression genes of ACE2 in colonic epithelial cells. Fig. 10B displayed that QFPD had some common targets with these five sets, while specific 254 targets for QFPD, indicating other mechanisms of QFPD on COVID-19 in addition to 2019-nCov, pneumonia, ACE2 and CD147 related functions.

Fig. 10

ACE2 and CD147 expression across tissues and co-expression genes. A: Radar plot of ACE2 and CD147 expression across 53 tissues. The expression values were converted to base-2 logarithm. Red triangle and square mean the top 5 expression tissues. B: UpSet plot of proteins among QFPD, HCoV (Host_protein), pneumonia, ACE2 co-expression genes (ACE2_database), CD147 co-expression genes (CD147_database), and ACE2 co-expression genes in colonic epithelial cells (ACE2_colonic). The horizontal bar graph at the bottom left shows the total number of proteins for each group set. Circles and vertical lines in the x-axis mark the corresponding data sets being compared. The vertical bar graph at the top quantitates the number of proteins in the comparisons. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Firstly, the robustness of whole networks against formula attack was assessed to evaluate QFPD attack on the COVID-19 disease network. Interestingly, Table 4 and Fig. 11 showed that MSXG, SGMH, XCH, WLS and Others attack on the COVID-19 network were characterized by greater disturbance score than negative control (BXTM), and increasing dependence on hub nodes, indicating greater fragility under formula attack. In addition, SGMH, MSXG, and Others exerted higher disturbance score than the positive control (YDBF) Table 5 .

Fig. 11

Evaluation of the effect of QFPD on the robustness disturbance of COVID-19 network. Blue normal distribution: drug attack on random networks as a null distribution for the permutation test. Red vertical line: the disturbance rate of the drug to the real disease network. Fist row: MXSG, second row: SGMH, third row: XCH, fourth row: WLS, fifth row: Others. Fist column: average connectivity, second row: average length of shortest path, third row: connection centrality, fourth row: closeness centrality. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 5

Docking score between common ingredients of QFPD and 2019-nCov proteins.

Molecule	MS1	MX16	MX17	MX8	MXO1	SO1	SX1	WO1	XO1
Main Protease	−7.6	−7	−7.8	−7.8	−7.8	−7.2	−7.1	−7.4	−7.4
Papain-like protease	−8	−8.2	−8.1	−8.1	−8.2	−8.2	−9.2	−8.7	−8.3
RdRp with RNA	−8.6	−8	−7.9	−7.9	−8	−8.4	−8.2	−8.8	−9.5
RdRp without RNA	−7.1	−6.8	−6.8	−6.8	−7	−7.1	−6.9	−7.5	−7.1
Helicase ADP site	−6.2	−6.5	−6.3	−6.1	−6.2	−6.4	−6.8	−6.2	−6.8
Helicase NCB site	−7.4	−7.1	−7.2	−7.2	−7.5	−7.5	−7.4	−7.6	−7.5
Nsp14(ExoN)	−7.2	−6.6	−6.6	−6.7	−6.9	−7	−6.9	−7.1	−7.1
Nsp14(N7-MTase)	−8.6	−8.3	−8.4	−8.3	−8.6	−8.3	−9.4	−8.7	−8.7
Nsp15(endoribonuclease)	−6.2	−6.8	−6.4	−6.4	−6.3	−6.5	−6.8	−6.4	−6.6
Nsp16(2′-O-MTase)	−8.1	−7.4	−7.6	−7.8	−7.8	−8.3	−7.7	−8.4	−8.2
N protein NCB site	−8.1	−7.9	−8	−8	−7.6	−7.8	−7.5	−7.7	−7.9
E protein(ion channel)	−6.4	−6.9	−7.2	−7.2	−6.9	−6.9	−7.9	−6.8	−6.8

M: MXSG, S: SGMH, X: XCH, O: Others.

Next, to illustrate the mechanism of QFPD against COVID-19, a formula-attacked target-KEGG pathway network was constructed (Fig. 12 ). This network showed that MSXG, SGMH, XCH, WLS and Others interacted with 8 drug-attacked nodes (Cdc20, Ido1, Ifng, Il10, Il6, Ptger4, Spi1, Tnf), and 24 drug-attacked nodes in the COVID-19 network were related to Graft-versus-host disease, cytokine-cytokine receptor interaction, asthma, influenza A, inflammatory bowel disease, JAK-STAT signaling, etc.

Fig. 12

Disturbance analysis of QFPD for COVID-19 network. A: Venn diagram of the five formulae’ attacked targets. B: Formula-attacked target-KEGG pathway network; green square: formula, yellow diamond: attacked target, red circle: KEGG pathways. The bigger the size of the nodes is, the higher the degree is. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Discussion

Novel coronavirus (2019-nCov) infection is characterized by lung and immune system damage. Severe infection can lead to acute respiratory distress syndrome (ARDS) and septicemia, and eventually lead to death [2]. In addition, a number of patients presented multi-organ damage and dysfunction [28]. However, there are no specific drugs or vaccine for the treatment of the COVID-19. This reason may partly be that a single targeted drug cannot cure a complex disease with complex biological networks [29]. Despite the lack of strong evidence-based medicine, TCM has a good potential to complement the medical service for COVID-19, including reverting radiological changes, and shortening fever duration and hospital stay. [30]. It is observed that the total effective rate of QFPD in the treatment of pneumonia patients infected by novel coronavirus is more than 90 % [31]. Therefore, we explored the mechanism of QFPD against COVID-19 by systems pharmacology, and provided a combination strategy to explore the functional units in QFPD from a holistic perspective.

To our knowledge, this is the first study to explore the mechanisms of QFPD for COVID based on intra-functional units. In this study, GO enrichment analysis showed that the common GO terms of MSXG, SGMH, XCH, WLS and Others targets were significant enriched in oxidoreductase activity, lipid metabolic process, lipid binding, small molecule metabolic process, and homeostatic process etc., suggesting QFPD may exert anti-viral activity through metabolic function. In agreement with these results, a recent research has found that lipid metabolic reprograming plays an important role in virus replication, which may be an appealing and applicable target for antiviral therapy [32]. KEGG analysis showed that in addition to lipid metabolism-related pathways, endocrine system pathways were also significantly enriched in more than four formulae, including PPAR signaling pathway and adipocytokine signaling pathway. A recent study has showed that the host can exert anti-inflammatory functions to inhibit excessive inflammatory damage through PPAR signaling pathway after H1N1 infection, thus keeping homeostasis of metabolism and development [33]. In addition, other common terms were significant enriched in more than two formulae, such as immune system process, endoplasmic reticulum, cell-cell signaling, calcium signaling pathway, vascular smooth muscle contraction, inflammatory mediator regulation of TRP channels, cardiac muscle contraction, etc. Therefore, the multi-pathway and multi-target results of our intra functional unit of QFPD not only showed a new useful method for studying TCM, but may demonstrate the rationality of TCM compatibility. Moreover, TCMATCOV platform was used to validate these results. Interestingly, all the five FUs of QFPD showed higher disturbance score than negative control (BXTM), indicating that MSXG, SGMH, XCH, WLS and Others may protect COVID-19 independently, and target 8 specifically expressed drug-attacked nodes (Cdc20, Ido1, Ifng, Il10, Il6, Ptger4, Spi1, Tnf) which were related to the bacterial and viral responses, cytokine, immune system, signaling transduction, etc.

Currently, a number of studies have showed that 2019-nCov can cause multiple organs dysfunction, including liver [34], pancreas [35], kidney [36], throat and rectum [37], which may be the reason that a wide distribution of ACE2 across these tissues [38]. In agreement with these results, we found that ACE2 was highly expressed in terminal ileum, testis, adipose visceral omentum, heart left ventricle, kidney cortex and thyroid, etc., and QFPD has only 15.33 % common targets with 2019-nCov, pneumonia and ACE2 related genes. These results indicate that the effective treatment of QFPD for COVID-19 may be through a holistic treatment. Moreover, TTD analysis further displayed that QFPD targets were significantly enriched in many COVID-19 related disease, such as chronic inflammatory diseases, asthma, inflammatory bowel disease, chronic obstructive pulmonary disease, intrahepatic cholestasis, chronic ileitis, etc.

It is known that ADMETox prediction is an important part in evaluating if a drug can be toxic or can be absorbed during drug development process [39]. In our study, ADMETox evaluation shows that 20 compounds passed the stringent lead-like criteria (250 ≤ MW≤350 & XLOGP ≤ 3.5 & Number of rotatable bonds≤7) [40] and in silico drug-likeness test, and showed high gastrointestinal absorption. Moreover, predicted toxicity evaluation showed that the median lethal dose (LD50) of all these ingredients was above 1600 mg/kg, thus may suggesting safety and efficacy of QFPD. Combined with molecular docking results, 4 specific ingredients (M3, S1, X2 and O2) and 5 common ingredients (MS1, MX16, SX1, WO1 and XO1) of QFPD might be promising leading compounds with good molecular docking score for 2019-nCov structure and non-structure proteins, revealing that QFPD treated COVID-19 by multi-component synergy. However, these newly monomer components should provide a further research.

It has been reported that host cellular microRNAs (miRNAs) are involved in the regulation of virus infection [41]. A previous study discovered that significantly up-regulated MIR301 and down-regulated MIR183/130B were found in H1N1 patients [42]. Consistent with these results, we found that MIR183 and MIR130A/B/301 are related to four functional units of QFPD, indicating these microRNAs may exert anti−COVID-19 activity through QFPD. In addition, CDKs have played a role in the efficient replication of various viruses, including human HIV-1, papillomaviruses, human cytomegalovirus (HCMV), herpes simplex virus (HSV) type 1 and HSV-2 [43,44]. In agreement with these results, we found that CDK7 was predicted to enriched in the five formulae, suggesting that QFPD may regulate replication of COVID-19 viruses via CDK7 mediated cell cycle and RNA polymerase II transcription. Recently, a previous study showed that LXR known to regulate cholesterol homeostasis during inflammation were differentially regulated during H1N1 influenza virus infection [45]. Based on our results that LXR was associated with MSXG, SGMH, XCH, WLS targets, we speculated that QFPD can regulate metabolic and pro-inflammatory processes to counter COVID-19 virus infection.

In summary, QFPD is effective in the treatment of COVID-19. However, some shortcomings in our study include lack of an in-depth study of predictive monomers and key targets and pathways, thus need further validation in vivo and in vitro. And the TCMATCOV platform uses SARS disease network, which is different from the COVID-19 disease network, and COVID-19-related cytokines are related to severe COVID-19 disease, so the results of platform analysis more reflect the potential efficacy of severe stage. Nevertheless, this study confirms that network pharmacology can help explore the mechanism of QFPD on the treatment of COVID-19 with time- and cost-saving. Moreover, based on our new FUNP analysis, we reveal that QFPD treat COVID-19 by a holistic treatment and multi-component synergy, and are further demonstrated by formula perturbation analysis. In addition, this study provides possible candidate monomers of QFPD and related miRNAs, kinases and TFs with potential therapeutic effect on COVID-19. This will hopefully provide evidence and new insights for further researches on the treatment of COVID-19 using QFPD.

Author contributions

Conceiving the research, Jian Chen and Yong-bing Cao; Data curation, Jian Chen, Wen-jie Sun, and Zhi-qiang Liang; Funding acquisition, Jian Chen, Zhi-qiang Liang, Bing-yong Cao and Ye-min Cao; Investigation, Ling-San Hu, Jian-ru Wang, and Bing-yong Cao; Methodology, Yong-kui Wang, Jiang-wei Yang and Ye-min Cao; Resources, Ling-San Hu and Yong-kui Wang; Visualization, Jian-ru Wang, Jiang-wei Yang; Writing – original draft, Jian Chen.

Declaration of Competing Interest

The authors declare no conflict of interest.