Mechanism of Dayuanyin in the treatment of coronavirus disease 2019 based on network pharmacology and molecular docking

Background At present, coronavirus disease 2019 (COVID-19), caused by infection with severe acute respiratory syndrome coronavirus 2, is spreading all over the world, with disastrous consequences for people of all countries. The traditional Chinese medicine prescription Dayuanyin (DYY), a classic prescription for the treatment of plague, has shown significant effects in the treatment of COVID-19. However, its specific mechanism of action has not yet been clarified. This study aims to explore the mechanism of action of DYY in the treatment of COVID-19 with the hope of providing a theoretical basis for its clinical application. Methods First, the TCMSP database was searched to screen the active ingredients and corresponding target genes of the DYY prescription and to further identify the core compounds in the active ingredient. Simultaneously, the Genecards database was searched to identify targets related to COVID-19. Then, the STRING database was applied to analyse protein–protein interaction, and Cytoscape software was used to draw a network diagram. The R language and DAVID database were used to analyse GO biological processes and KEGG pathway enrichment. Second, AutoDock Vina and other software were used for molecular docking of core targets and core compounds. Finally, before and after application of DYY, the core target gene IL6 of COVID-19 patients was detected by ELISA to validate the clinical effects. Results First, 174 compounds, 7053 target genes of DYY and 251 genes related to COVID-19 were selected, among which there were 45 target genes of DYY associated with treatment of COVID-19. This study demonstrated that the use of DYY in the treatment of COVID-19 involved a variety of biological processes, and DYY acted on key targets such as IL6, ILIB, and CCL2 through signaling pathways such as the IL-17 signaling pathway, AGE-RAGE signaling pathway in diabetic complications, and cytokine–cytokine receptor interaction. DYY might play a vital role in treating COVID-19 by suppressing the inflammatory storm and regulating immune function. Second, the molecular docking results showed that there was a certain affinity between the core compounds (kaempferol, quercetin, 7-Methoxy-2-methyl isoflavone, naringenin, formononetin) and core target genes (IL6, IL1B, CCL2). Finally, clinical studies showed that the level of IL6 was elevated in COVID-19 patients, and DYY can reduce its levels. Conclusions DYY may treat COVID-19 through multiple targets, multiple channels, and multiple pathways and is worthy of clinical application and promotion.


Background
Coronavirus disease 2019 (COVID- 19) was first discovered in Wuhan, China, on December 12, 2019, but to date, no definitive conclusion has been drawn about its origin. According to the classification of syndromes in traditional Chinese medicine, COVID-19 is classified as an "epidemic disease" (damp-warm disease), and it is a highly contagious disease. In the early stages of dampwarm diseases, "damp-warm disease with syndrome of pathogen blocking pleuro-diaphragmatic interspace" is very common and is a specific stage and phenomenon in the pathological process of the disease. This symptom first appeared in the Theory of Epidemic Febrile Disease by Wu Youke during the Ming Dynasty, and he created Dayuanyin (DYY), described in the book. After 2019-nCoV invades the human body, it disturbs and damages the human immune system, further causing different degrees of damage to various organs throughout the body [1,2]. DYY, a traditional Chinese medicine prescription, has played an important role in the prevention and treatment of epidemic diseases in documented history and literature. It has been used to treat influenza [3], atypical pneumonia [4], AIDS [5] and other diseases and has proven to be very effective in clinical applications. At the same time, through clinical observation of COVID-19 patients in the early stage of DYY treatment, it was found that this prescription can improve the clinical symptoms and signs of patients, improve the prognosis of patients, and shorten the course of disease [6,7], making it worthy of clinical application and promotion. However, its mechanism of action in COVID-19 patients has not yet been clarified.
Network pharmacology, originally proposed by Andrew L Hopkins, includes systems biology, pharmacology, mathematics, computer network analysis, etc. As a useful tool for systematically evaluating and demonstrating the rationality of drugs, it has now been widely accepted [8,9]. The application of network pharmacology in traditional Chinese medicine provides us with new possibilities for screening active ingredients of drugs and targets for disease treatment, which is helpful for explaining the mechanism of action of drugs for disease treatment at a system level [10]. Molecular docking is a theoretical simulation method that mainly studies intermolecular interactions and predicts their binding mode and affinity [11]. Not only can it be used for drug development, but it can also provide keen insights into protein function prediction and other important issues [12].
This study aimed to use network pharmacology and molecular docking to preliminarily explore the mechanism of action of this prescription in the treatment of COVID-19 patients, with the goal of widely using this prescription for COVID-19 patients with early dampwarm syndromes to improve the patients' condition and to prevent the ongoing COVID-19 outbreak. A technological road-map of the experimental procedures of our study is shown in Fig. 1.

Acquisition of the chemical composition and target information of DYY
The Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) records 499 common traditional Chinese medicines (Chinese Pharmacopoeia 2010 edition) and elaborates their ingredients, the corresponding target information and common disease information related to traditional Chinese medicines [13]. The database provides pharmacokinetic information for each compound, such as drug-like (DL), oral bioavailability (OB), and blood-brain barrier (BBB). In this study, the TCMSP database (http://tcmsp w.com/ tcmsp .php) was used to search for and determine the active ingredients in the composition of the DYY decoction. At the same time, target genes were predicted for these active ingredients. OB and DL property are important reference standards for evaluating whether compounds can be used as drugs. In this study, OB ≥ 30% and DL ≥ 0.18 were used as screening thresholds [14]. According to the selected active ingredients of DYY, the target genes corresponding to the above active ingredients derived from DrugBank were further screened using Perl language in combination with the TCMSP database.

Construction of network diagrams of compounds and corresponding targets
The compounds and predicted targets in the DYY formula obtained through the TCMSP database were imported into Cytoscape 3.6.1 software, and a compound-target network diagram was drawn to obtain the top five core compounds.

Acquisition of disease targets
The keyword "novel coronavirus pneumonia" was entered into the Genecards (https ://www.genec ards.org/, version 4.12) database to obtain target genes related to the COVID-19 disease.

Intersection of disease genes and drug genes
The target genes predicted from the active ingredients in DYY were intersected and mapped with the target genes predicted for the COVID-19 disease to obtain the target genes of DYY for the treatment of COVID-19. The Venn Diagram package in R was used to draw a Venn diagram.

Protein-protein interaction analysis and core target screening
The target genes for DYY treatment of COVID-19 were entered into the STRING database for protein-protein interaction (PPI) analysis, "Homo sapiens" was selected, the minimum required interaction score was set to > 0.9, the protein interaction network map was downloaded, and R3.5.0 was used to screen core genes.

Construction of network visualization
The active ingredients of DYY, the targets corresponding to the active ingredients, and the targets predicted for the COVID-19 disease were imported into Cytoscape 3.6.1 software, and a drug-target-disease network diagram was constructed for network visualization.

GO analysis and KEGG pathway enrichment analysis
The DAVID database (http://david .abcc.ncifc rf.gov/) can functionally annotate many genes and help us understand the biological process and meaning behind genes. The target genes selected above were combined with R language and DAVID database for Gene Ontology (GO) biological process enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment.

Construction of KEGG relationship network
The pathway ID numbers and the genes involved in the KEGG-enriched pathways were imported into Cytoscape software, 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 to construct a KEGG relationship network.

Molecular docking verification of core compounds and core target genes
Firstly, the top five core compounds were selected, and the two-dimensional structure diagrams of the compounds were downloaded from the PubChem database, imported into Chem3D software to draw the threedimensional structure diagrams of the core compounds and optimize the energy, and saved in mol2 format. Then, the files were imported into AutoDockTools-1.5.6 software to add charge and display rotatable keys and then saved in pdbqt format. Secondly, the protein crystal structures corresponding to the core target genes were downloaded from the PDB database, imported into Pymol software to remove water molecules and heteromolecules, imported into AutoDockTools-1.5.6 software to add hydrogen atoms and charge operations, saved to pdbqt format, and imported into Discovery Studio 3.5 Client software to search for active pockets. Finally, the above core compounds were used as ligands, and the proteins corresponding to the core target genes were used as receptors for molecular docking. The results were analysed and interpreted using PyMOL software and Discovery Studio 3.5 Client.

Clinical validation of the core target IL6
In this study, a total of 45 patients who were hospitalized in Third People's Hospital of Hubei Province and Lei Shen Shan Hospital during the period from January 25, 2019 to March 8, 2020 were selected. The TCM syndrome of the selected patients was "plague (syndrome of pathogen hidden in interpleuro-diaphramatic space)", and DYY was used for treatment. ELISA was used to detect changes in IL6 levels at the time of before and after treatment for 1 week. All statistical analyses were performed by GraphPad Prism version 7.00 software. T test was used for comparisons before and after treatment in each group. Our data are expressed as the mean ± standard deviation (SD). A value of p < 0.05 was considered significant.

Mechanism of action of DYY in the treatment of COVID-19
Adobe Illustrator CC software was used to draw the chart for the specific mechanism of DYY treatment of COVID-19.
According to the results obtained by screening the active ingredients against the TCMSP database, there were a total of 7053 targets in the DrugBank. Among them, the number of targets of magnolia officinalis, amomum, arecae semen, herbaceous peony, scutellariae radix, anemarrhenae rhizoma and licorice was

Construction of network diagrams of compounds and corresponding targets
The compounds and corresponding targets in the DYY formula were imported into Cytoscape software to draw a compound-target network diagram (see Fig. 2). In this study, degree was selected as a measure of node importance. With the help of the Network Analyzer plug-in in Cytoscape software, the topology parameters of the network were calculated and analysed from the perspective of network node importance. Degree refers to the number of edges associated with a node. The greater the degree of a node is, the larger the node area in the graph. That is, the larger the node area is, the greater the importance of the node in the network. The compounds in Fig. 2 and their corresponding targets were used as network nodes. Figure 2 shows that one compound can act on multiple target genes, and multiple compounds can also act on one target gene at the same time. Among the compounds, MOL000422 (kaempferol), MOL000098 (quercetin), MOL003896 (7-Methoxy-2-methyl isoflavone), MOL004328 (naringenin), MOL000392 (formononetin) and MOL000358 (beta-sitosterol) occupied the largest area on the graph among all compounds and were important core compounds.

Acquisition of disease target
A total of 251 genes related to COVID-19 were obtained by searching the Genecards database. The relevance score was used as the selection criterion to obtain the top 30 genes (see Table 2).

Intersection of drug targets and disease targets
The above DYY drug target genes were intersected with COVID-19 disease targets to obtain possible genes associated with DYY treatment of COVID-19. The results showed that there was a total of 45 genes associated with DYY treatment of COVID-19 (see Additional file 2: Fig.  S2).

PPI analysis and core target screening
The STRING database was used to draw a PPI network diagram of DYY for COVID-19 (see Fig. 3a). As shown in Fig. 3a, the network diagram consisted of 45 nodes and 581 edges, for which the average node degree was 25.8, and the PPI enrichment p-value was < 1.0e−16. The above PPI network was processed using R language, and the top 30 core genes were selected (see Fig. 3b). Figure 3b shows that the top 30 core genes had a node degree greater than 21, and the top genes, such as IL6, MAPK3, MAPK8, CASP3, IL10, IL1B, CXCL8, MAPK1, CCL2, IFNG and IL4, had a higher number of connections than other genes, all showing 35 or more connections.

Construction of network visualization
The active ingredients of DYY, the targets corresponding to the active ingredients, and the targets predicted for COVID-19 were imported into Cytoscape software to build a drug-target-disease network diagram (see Fig. 4). The network had a total of 139 nodes (including 94 compound nodes and 45 gene nodes) and 546 connections.

GO analysis and KEGG pathway enrichment analysis
The R language and DAVID database were used for GO enrichment analysis by using the above-mentioned targets of DYY to treat COVID-19, and the number of biological process (BP), cellular component (CC), and molecular function (MF) entries was 1,506, 33 and 83, respectively. The top 30 biological processes were screened and are represented as graphical bubbles (see Fig. 5a-c). The KEGG pathway enrichment analysis identified 40 signaling pathways, and the top 20 entries were selected and are represented by a bar graph (see Fig. 5d).

Construction of the KEGG relationship network
The top 20 pathways involved in DYY treatment of COVID-19 and the genes enriched in these pathways were imported into Cytoscape software to build a KEGG relationship network diagram (see Fig. 6). In Fig. 6, the pathways and the genes enriched in the pathways were used as network nodes, and we selected the top three pathways (hsa04657 (IL-17 signaling pathway), hsa04933 (AGE-RAGE signaling pathway in diabetic complications), and hsa0406 (cytokine-cytokine recev] ptor interaction pathway)) and top three target genes (IL6, ILIB, CCL2) enriched in these pathways according to degree.

Molecular docking verification of core compounds and core target genes
The results obtained by the molecular docking software are shown in Table 3. The letters x, y and z were used to represent the size and position of the pocket. The final selected pocket is shown in bold in the column titled 'Pocket size' . The results of the docking of the receptor and ligand are shown under 'Docked complex' , and the residues docked with the small-molecule ligand are shown as yellow sticks. The structure with the initial ligand and the predicted protein pocket were processed by Discovery Studio 3.5 Client software, and the docked complex was processed by PyMOL software. As seen from Table 3, the scores for the five core compounds (kaempferol, quercetin, 7-Methoxy-2-methyl isoflavone, naringenin, formononetin) and protein crystal structures corresponding to the core target genes (IL6, IL1B, CCL2) were all greater than -5 kcal/mol, indicating that the compound had a certain affinity for the protein crystal structure. The interactions between some    Fig. S3 shows that the small-molecule compounds were tightly bound to the protein residues via various interactions.

Clinical validation of the core target IL6
Of the 45 patients selected, 15 were mild cases (group A), 15 were moderate cases (group B), and 15 were severe cases (group C). The age distribution of each group of patients and changes in IL6 levels before and after treatment are shown in Fig. 7a, b, respectively. Compared with the severe cases, Fig. 7a shows that the mild and moderate cases were younger. Figure 7b shows that a majority of the patients had different levels of IL6 elevation before treatment (the normal reference value of IL6 is 0-7 pg/ ml), and the increase in IL6 was most pronounced in severe cases. After treatment, IL6 decreased in all groups, and differences within each group before and after treatment were statistically significant.

Mechanism of action of DYY in the treatment of COVID-19
Based on the above studies, the specific mechanism of the action of DYY in the treatment of COVID-19 is shown in Fig. 8.

Discussion
As of April 6, 2020, the cumulative number of COVID-19 confirmed cases worldwide has exceeded 1.2 million. However, no vaccine or definitive antiviral drugs are available for the prevention and treatment of COVID-19. Therefore, it is crucial to find out medicines with confirmed curative effects for the prevention and treatment of COVID-19 as soon as possible to improve the patient's condition and prevent the ongoing outbreak of COVID-19.
In this study, we first searched and screened a database of traditional Chinese medicine database to obtain 174 DYY compounds and 7053 corresponding target genes. Ent-epicatechin, quercetin, mairin, beta-sitosterol, sitosterol, and stigmasterol are common compounds of two or more Chinese medicines. Studies have shown that quercetin can reduce apoptosis induced by hypoxia and rescue phosphorylation of AMPK [15]. Beta-sitosterol has antipyretic, analgesic, anti-inflammatory, and antioxidant functions and plays roles in cough and phlegm elimination, immune regulation and tissue repair [16,17]. Stigmasterol is present in the membrane [18] and has anti-inflammatory effects [19]. The compound-target network diagram (Fig. 2) shows that there was a complex network relationship between the compounds and the targets. Kaempferol, quercetin, 7-Methoxy-2-methyl isoflavone, naringenin, formononetin, and beta-sitosterol had the largest number of targets and were core compounds.
By searching the disease database, we found a total of 251 genes related to COVID-19. Drug target genes and disease related genes were intersected, and a total of 45 target genes for DYY treatment of COVID-19 were finally obtained (Additional file 2: Fig. S2). These genes were analysed by PPI analysis to obtain the corresponding network diagram (Fig. 3a). Figure 3a shows that the target genes of DYY for the treatment of COVID-19 were not independent, and there was a certain relationship among these genes. The core gene map in Fig. 3b shows that the top-ranking genes were IL6, MAPK3, MAPK8, CASP3, IL10, IL1B, CXCL8, MAPK1, CCL2, IFNG, IL4, etc. These genes were mainly concentrated in the inflammatory response, immune modulation, and cellular stress processes, which indicated that they might play a key role in DYY treatment of COVID-19. It is well known that IL6, IL10, IL1B, and IL4 are all members of the interleukin family. Interleukins play an important role in transmitting information, regulating immune cells, mediating T and B cell activation, and responding to inflammation [20,21]. CCL2 and CXCL8 belong to the chemokine family and are important inflammatory cytokines. They play an important role in the migration of Tregs to inflammatory tissues [22] and in immune regulation in the body  (COVID-19), the green arrows represent the compounds, and the pink boxes represent the targets of action [23]. MAPK3, MAPK8, and MAPK1 are members of the MAPK family and can participate in responses to potentially harmful abiotic stress stimuli [24].
At the same time, a network of drug active ingredients, target genes corresponding to the active ingredients and disease targets was constructed as shown in Fig. 4. We can see that one compound can act on multiple target genes. Similarly, one target gene can also correspond to multiple compounds. That is, multiple compounds can act on a common target. Based on the above analysis, we concluded that multiple active ingredients in the traditional Chinese medicine prescription DYY can act on COVID-19 through multiple targets.
Through functional enrichment analysis of target genes for DYY treatment of COVID-19, GO biological process and KEGG pathway enrichment maps were obtained (see Fig. 5). It can be seen from Fig. 5 that in the GO terms, the BP terms (Fig. 5a) were mainly associated with the cell's response to processes such as lipopolysaccharide, molecule of bacterial origin, biotic In the bubble charts in (a-c), the ordinate represents the names of the BP, CC, and MF terms, respectively, and the abscissa represents the degree of enrichment. In (d), the ordinate represents the names of the pathways, and the abscissa represents the number of genes enriched in the pathway. The P value indicates the significance of enrichment. The smaller the P value is, the higher the significance of enrichment, and the redder the colour on the graph  stimulus, cytokine production, oxidative stress, and adaptive signaling pathways. CC terms (Fig. 5b) were mainly associated with various membranes, including membrane raft, membrane microdomain, membrane region, plasma membrane raft, nuclear envelope, and mitochondrial outer membrane, etc.; the terms were also associated with focal adhesion, cell-substrate adherens junction, cell-substrate junction, etc. The MF terms (Fig. 5c) were mainly associated with various receptors (cytokine, chemokine, growth factor, CXCR chemokine, G protein-coupled, nuclear hormone receptors), binding functions (phosphatase, heme, protein phosphatase, tetrapyrrole, BH domain, death domain, tumor necrosis factor receptor superfamily) and various cytokine, ligand, and kinase activities (cytokine, receptor ligand, MAP kinase, chemokine, oxidoreductase). The pathways involved in the KEGG enrichment pathway (Fig. 5d) Fig. 6, hsa04657 (IL-17 signaling pathway), hsa04933 (AGE-RAGE signaling pathway in diabetic complications), and hsa0406 (cytokine-cytokine receptor interaction pathway) enriched the highest number of genes, indicating that these pathways may play an important role in the mechanism of action of DYY in the treatment of COVID-19. The IL-17 signaling pathway is involved in the body's immune response [25,26] and inflammatory response [27]. The AGE-RAGE signaling pathway has important protective effects on bones and the heart and participates in oxidative stress response [25] and fibrosis transduction [22]. The cytokine-cytokine receptor interaction is a key pathway for regulating the cellular inflammatory response [28]. IL6, ILIB, CCL2 and other genes occupy a large rectangular area, indicating that there are more pathways connected to these genes. Therefore, it can be speculated that these genes play a key role in the mechanism of action of DYY in the treatment of COVID-19. IL6, ILIB, and CCL2 represent a wide range of inflammatory mediators and pathways. Many animal and human experiments have demonstrated that IL6 has a wide range of anti-inflammatory effects [29]. ILIB has analgesic, immunomodulatory, anti-hypoxia, and anti-inflammatory functions. CCL2 is an important inflammasome-associated chemokine [30]. Inhibition of CCL2 can reduce the infiltration of peripheral inflammatory cells such as monocytes and neutrophils [9]. NOS3 is a vasoprotective gene [31] that regulates vascular tone, blood pressure and platelet aggregation [32]. Research reports have shown that NOS3 can affect metabolism in the urea cycle of the methylation pathway, which is essential for preventing systemic inflammation [33]. Graph of patient age distribution and IL 6 levels. In (a), the abscissa represents group, and the ordinate represents age. In (b), the abscissa represents the group, and the ordinate represents the level of IL6, where A1, B1, and C1 represent patients before treatment while A2, B2, and C2 represent patients after treatment in each group. Data are expressed as the mean ± S.D. (n = 15 per group); *p < 0.01, **p < 0.001, and***p < 0.001. the after treatment groups(A2, B2, C2) vs. the before treatment groups(A1, B1, C1) By combining the core target gene bar chart (Fig. 3b) and the KEGG relationship network diagram (Fig. 6), we can see that IL6 is one of the most critical genes for anti-inflammatory and immune regulation in COVID-19 patients treated with DYY. Based on the comparison of COVID-19 patients before and after treatment with DYY, the IL6 level of COVID-19 patients increased to different degrees when they were admitted to the hospital but decreased after treatment, further confirming that DYY may play an important role in anti-inflammatory and immune regulation and may have other effects in the treatment of COVID-19 patients.

Conclusions
In summary, we speculate that DYY may play an antiinflammatory and immunoregulatory role in COVID-19 by acting on multiple target proteins, such as IL6, ILIB, and CCL2. The role of DYY involves a variety of biological processes, mainly signaling pathways such as the IL-17 signaling pathway, cytokine-cytokine receptor interaction, and AGE-RAGE signaling pathway, involved in diabetic complications. In short, DYY