Exploration of molecular targets and mechanisms of Chinese medicinal formula Acacia Catechu -Scutellariae Radix in the treatment of COVID-19 by a systems pharmacology strategy.

Coronavirus disease 2019 (COVID-19) is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). In China, the Acacia catechu (AC)-Scutellariae Radix (SR) formula has been widely used for pulmonary infection in clinical practice for several centuries. However, the potential role and mechanisms of this formula against COVID-19 remains unclear. The present study was designed to dissect the active ingredients, molecular targets, and the therapeutic mechanisms of AC-SR formula in the treatment of COVID-19 based on a systems pharmacology strategy integrated by ADME screening, target prediction, network analysis, GO and KEGG enrichment analysis, molecular docking, and molecular dynamic (MD) simulations. Finally, Quercetin, Fisetin(1-), kaempferol, Wogonin, Beta-sitosterol, Baicalein, Skullcapflavone II, Stigmasterol were primarily screened to be the potentially effective active ingredients against COVID-19. The hub-proteins were TP53, JUN, ESR1, MAPK1, Akt1, HSP90AA1, TNF, IL-6, SRC, and RELA. The potential mechanisms of AC-SR formula in the treatment of COVID-19 were the TNF signaling pathway, PI3K-Akt signaling pathway and IL-17 signaling pathway, etc. Furthermore, virtual docking revealed that baicalein, (+)-catechin and fisetin(1-) exhibited high affinity to SARS-CoV-2 3CLpro, which has validated by the FRET-based enzymatic inhibitory assays with the IC50 of 11.3, 23.8, and 44.1 μM, respectively. And also, a concentration-dependent inhibition of baicalein, quercetin and (+)-catechin against SARS-CoV-2 ACE2 was observed with the IC50 of 138.2, 141.3, and 348.4 μM, respectively. These findings suggested AC-SR formula exerted therapeutic effects involving "multi-compounds and multi-targets." It might be working through directly inhibiting the virus, improving immune function, and reducing the inflammatory in response to anti-COVID-19. Ultimately, this study would provide new perspective for discovering potential drugs and mechanisms against COVID-19.

INTRODUCTION

The number of confirmed infections has climbed to more than 458.5 million and approximately 6.1 million people dead globally since the outbreak of corona virus disease 2019 (COVID‐19) in December 2019 (https://covid19.who.int/), caused by the severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) (Lipsitch, Swerdlow, & Finelli, 2020). Infection with the SARS‐CoV‐2 may result in lung inflammation and infiltration, with the frequent manifestations including fever, coughing, and shortness of breath (Chen et al., 2020; Merza et al., 2021). Severe instances could develop into acute respiratory distress syndrome (ARDS), which results in systemic inflammatory cytokine storms, with unpredictably negative effects for the outcome of the disease (Richardson et al., 2020; Wiersinga, Rhodes, Cheng, Peacock, & Prescott, 2020). The worldwide spread of COVID‐19 has posed a profound threat to human health, and the implementation of vaccines is still a major asset in slowing down the pandemic by far (Fiolet, Kherabi, CJ, Ghosn, & Peiffer‐Smadja, 2022; Huang, Bai, He, Xie, & Zhou, 2020). However, vaccines are slightly less effective against worrisome variants like as Delta and Omicron, while there are remain few specific drugs to cure SARS‐CoV‐2 infection on the market as of the end of 2021 (Brüssow, 2021; Mohammadi, Shayestehpour, & Mirzaei, 2021). Therefore, the global hunts for potential pharmaceuticals in adequately managing this disease are extremely urgent. In this regard, medicinal plants containing specific phytomoieties may provide a broad range of potential therapeutic applications (Anand et al., 2021; Brendler et al., 2021; Das et al., 2021).

Traditional Chinese Medicine has a history of almost 3,000 years in the prevention and treatment of infectious diseases. It has demonstrated extraordinary benefits in preventing and treating viral respiratory diseases such as MERS, SARS (Hsu et al., 2006; Lau et al., 2005), and H1N1 (Wang et al., 2011) in recent years and provides a number of distinct features, including preventive treatment of disease, therapy based on syndrome differentiation, and multitarget intervention (Xian et al., 2020). During the fight against the COVID‐19, TCM has also made great contributions in improving cure rate, shortening illness duration, delaying disease progression, and reducing mortality rate (Ren, Zhang, & Wang, 2020; Shahrajabian, Sun, Soleymani, & Cheng, 2021; Zhang et al., 2021). According to syndrome differentiation‐based therapy theory of TCM, the pandemic falls under the classification of “damp epidemic”. The damp epidemic pathogen prefers and is mostly located in the lung, with the sickness affecting the spleen and stomach as well as the liver and kidneys in extreme cases (Shi et al., 2020; Yuan, Xin, Tang, & Cong, 2020). The pathophysiology is centered on a dampness‐toxin obstructing the lung and suppressing Qi, which manifests as dampness, heat, poison, blood stasis, and deficiency (Qiu et al., 2020; Ren et al., 2021). As a result, Chinese physicians regard “expelling evil and detoxifying, drying moisture, and removing blood stasis as the most significant aspect” in the face of COVID‐19.

Acacia catechu (AC), also called Er‐Cha in Chinese, is a traditional medicinal plant having antitussive, antipyretic, hemostatic, and hepatoprotective properties that is widely used in China, India, and Southeast Asia (Khare, 2007). Scutellariae Radix (SR), also known as Huang‐Qin in Chinese, is a traditional Chinese herbal medicine that has historically been used to treat respiratory inflammatory or viral infections which cause throat swelling and soreness, pneumonia, and fever (Song et al., 2020). Actually, in clinical practice over the last several centuries, AC and SR have been regularly recommended in combination, famous as Huang‐qin Er‐cha Decoction, for the prevention or treatment of cough, phlegm, and fever caused by pulmonary infection in China (Wang et al., 2019). The combinational prescript of the two medications was believed to enhance their preventative or therapeutic benefits in the treatment of pulmonary disease. As demonstrated in our previous studies, AC‐SR formula exerted a strong anti‐inflammatory effect in LPS‐induced ALI. AC‐SR formula dramatically lowered the wet‐to‐dry weight ratio of the lungs, ameliorated LPS‐induced lung histopathological alterations, decreased inflammation, blunted the production of proinflammatory cytokines, such as interleukin‐1 (IL‐1) and tumor necrosis factor‐α (TNF‐α). The putative molecular mechanism for the protective effect of AC‐SR formula against ALI is responsible for the attenuation of the NF‐B, MAPKs, and PI3K‐Akt signaling pathways in LPS‐induced AEC‐II (Feng et al., 2019). In addition, according to the literature studies, SR as well as flavonoids from AC, such as kaempferol, quercetin, and (+)‐Catechin, could possesses potential antiviral and anti‐inflammatory effect against COVID‐19 (Jena, Kanungo, Nayak, Chainy, & Dandapat, 2021; Khazdair, Anaeigoudari, & Agbor, 2021; Liu et al., 2021; Song et al., 2020). Currently, however, there is little study focusing on the utilization of the AC‐SR formula in the treatment of COVID‐19. As we have reason to speculate that the anti‐SARS‐CoV‐2 mechanism of AC‐SR formula is related to improving body immunity and resisting cytokine storm based on the characteristics of this prescription, it is thus critical to discover an appropriate approach to explain the precise mechanism of this formula.

The application of systems biology methods to the research of SARS‐CoV‐2 is not only beneficial in deciphering the pathogenesis, as well as the molecular interactions that occur during infection, but also help to develop novel treatment strategies for the COVID‐19 pandemic (Banaganapalli et al., 2021; Wynants et al., 2021). Network pharmacology has been offered as a promising method to dissect herbal formulas and predict potential new drugs or targets for the COVID‐19 (Hong, Duan, Wu, Yang, & Wu, 2020; Tao et al., 2021; Xia et al., 2020), while molecular docking as well as molecular dynamics simulations represent an unique avenue for structural molecular biology and the computer aided drug design in the development of novel medications (Aljarba, Hasnain, Bin‐Meferij, & Alkahtani, 2022; Morris & Lim‐Wilby, 2008; Saikia & Bordoloi, 2019; Wang et al., 2021). In this study, a systems pharmacology strategy (Figure 1) integrated by network pharmacology, molecular docking, and molecular dynamic simulation were adopted to investigate the mechanism of action underlying the effectiveness of AC‐SR formula in COVID‐19 therapy. It would provide new perspective for discovering potential drugs and mechanisms against COVID‐19 and the developed strategy could also be able to serve as role models for the research and development of other natural medicines.

FIGURE 1

Flow chart of the present systems pharmacology strategy

MATERIALS AND METHODS

Identification and screening of active compounds

Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP, https://tcmsp-e.com/), Traditional Chinese Medicines Integrated Database (TCMID, http://www.megabionet.org/tcmid/), and Bioinformatics Analysis Tool for Molecular mechANism of Traditional Chinese Medicine (BATMAN‐TCM, http://bionet.ncpsb.org/batman-tcm/) were used to authenticate all compounds of the AC‐SR formula. We selected compounds of AC‐SR formula according to the criterion of oral bioavailability (OB) ≥30% and drug‐likeness (DL) ≥ 0.18, which are the most important indicators for evaluating the characteristics of absorption, distribution, metabolism, and excretion (ADME) (Xu et al., 2012). In addition, we also search a large‐scale text and selected oral absorbable compounds with pharmacological activity, in order to supplement the compounds.

Identification of protein targets

The protein targets associated with active compounds were retrieved from the TCMSP database (https://old.tcmsp-e.com/index.php), which provided information of 6,511 drug molecules and 3,987 targets as well as the interactions between them (Ru et al., 2014). Then, Universal Protein Knowledgebase (Rolf et al., 2004) (UniProt, http://www.uniprot.org), an authoritative database of protein sequences, which comprised 54,247,468 sequence items was used to extract the targets, including the gene names and gene ID.

Predicting the targets of COVID‐19

We collected different genes associated with COVID‐19 from four resources (Fan et al., 2021). (1) Human Gene Database (GeneCards, (https://www.genecards.org/), (2) Therapeutic Target Database (TTD, http://db.idrblab.net/ttd/), (3) Comparative Toxicogenomics Database (CTD, http://ctdbase.org/), and (4) DisGeNET database (https://www.disgenet.org/). The search keywords were “Coronavirus Disease 2019,” “SARS‐COV‐2,” and “MERS‐COV”. The results of the above databases were integrated to obtain the target genes of COVID‐19 after deleting the repeated genes. To acquire candidate targets of AC‐SR formula acting on COVID‐19, we integrated the compounds' predicted targets of AC‐SR formula with target genes of COVID‐19 and chose those replicate genes.

Construction of protein–protein interactions (PPI) network

In order to further analyze the interaction between target proteins, the targets of AC‐SR formula and COVID‐19 were imported into the STRING database (https://string-db.org/) to build the PPI network interaction (Consortium UP, 2021). Cytoscape V3.7.1 was used to construct and visualize the PPI network (Shannon et al., 2003). “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

To investigate the functional annotation and involved pathways of genes. Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were implemented using Org.Hs.eg.db (Version 3.8.2) and ClusterProfiler (Version 3.9) packages in R (Version 3.6.2) (Tao et al., 2021). An adjusted p‐value of ≤ 0.05 was considered to identify the enricAC‐SR formula terms.

Component‐target molecular docking

The three‐dimensional (3D) structure of the key target protein was downloaded from the RCSB PDB protein structure database (https://www.rcsb.org/). AutoDock Vina (version 1.5.6) was used to remove the water molecules, isolate proteins, add the nonpolar hydrogen, and calculate Gasteiger charges for the structure (Du et al., 2022; Trott & Olson, 2010). Then, the corresponding format components and target proteins were signed in http://www.swissdock.ch/docking website to conduct the molecular docking experiment, and the docking binding energy [Estimated ΔG (kcal mol−1)] was predicted. According to the score of binding energy, the binding ability of active ingredients to certain target proteins was verified. Meanwhile, ritonavir, Nirmatrelvir, a potent SARS‐CoV‐2 3CLpro protease inhibitor (Reina & Iglesias, 2020), and SSAA09E2, an inhibitor of ACE2 (Bibi, Gul, Ali, & Kamal, 2021), were taken as positive control. The binding energy ≤0 kcal/mol indicated that the compound could bind and interact with the target, whereas the binding energy <−5 kcal/mol demonstrated a very strong binding force (Fan et al., 2021).

Molecular dynamic (MD) simulation

The molecular dynamics simulation study was performed using the Discovery Studio (DS) 2019 software package to evaluate the stability and interaction of ACE2 and 3CL receptors with Quercetin and Baicalein. The ligand‐receptor complex was placed in an orthorhombic box and solvated using an explicit periodic boundary solvated water model. Then sodium chloride was added to the system with the ionic strength of 0.145 to simulate the physiologic environment. Afterward, we subjected the system to the CHARMM force field (Li et al., 2021) and relaxed it through energy minimization (500 steps of conjugate gradient and 500 steps of steepest descent). And the final root means square gradient was 0.991. Then, the system’s temperature was slowly driven from an initial temperature of 50 K to the aimed temperature of 300 K within 4 ps. The time of equilibration simulations was 20 ps. Molecular dynamics simulation (production module) lasted for 200 ns with 1 ns time step. We completed the simulations under the normal pressure and the relatively constant temperature of nearly 300 K throughout the procedure. The particle mesh Ewald algorithm was applied for the calculation of long‐range electrostatics. And the linear constraint solver algorithm was adapted to identify all bonds involving hydrogen. The MD trajectory was written with 1,000 frames during the entire simulation run but only initial protein backbone frames were aligned to understand the stability of ligand‐protein complex. The root means square deviations (RMSD) and root mean square fluctuation (RMSF) were used to understand the stability of complex, which were essential to infer good binding affinities (Doniach & Eastman, 1999; Dubey, Tiwari, & Ojha, 2013).

3CLpro and ACE2 enzyme activity inhibition test

To further identify the possible inhibitory activities of the screened molecules against SARS‐CoV‐2 3CLpro, and ACE2, the compounds were tested by two commercial kits (enhanced 2019‐nCoV Mpro/3CLpro inhibitor screening kit, P0315S, Beyotime Biotech, Shanghai, China; ACE2 inhibitor screening kit, P0320S, Beyotime Biotech, Shanghai, China) based on fluorescence resonance energy transfer (FRET) protease assay (Xu et al., 2021), respectively. According to the operating instructions, ebselen and MLN‐4760 were chosen as the positive control, respectively. Enzyme activities were measured with saturated substrate concentration and different inhibitor concentrations. The enhanced fluorescence emission upon substrate cleavage was monitored at the excitation and emission wave lengths of 325 and 393 nm. The 50% inhibitory concentrations (IC50) were determined by plotting curves of percent inhibition versus compound concentration. Results are reported as IC50 values.

RESULTS

Active compounds of AC‐SR formula

In total, 126 compounds of AC‐SR formula were obtained from the TCMSP, TCMID, BATMAN‐TCM database, and relevant documents, including 33 from AC, 93 from SR. Besides, 89 active compounds were screened according to OB ≥30% and DL ≥0.18, including 23 from AC and 66 from SR. After eliminating repeated active compounds (MOL002914, MOL000073), 42 active compounds were retrieved. Information of some active compounds in AC‐SR formula is outlined in Table 1.

TABLE 1

Basic information of the potential effective ingredients of Scutellaria Radix and Acacia Catechu formula

No.	Mol ID	Chemical component	OB (%)	DL	Herb	Structure
1	MOL008420	5‐Hydroxy‐2‐[2‐(4‐hydroxyphenyl)acetyl]‐3‐methoxy‐benzoic acid	93.33832	0.20887	Acacia catechu
2	MOL000492	(+)‐catechin	54.82643	0.24164	Acacia catechu
3	MOL008428	3,4,8,10‐tetrahydroxy‐5H‐chromeno[3,2‐c]isochromen‐7‐one	50.53474	0.49886	Acacia catechu
4	MOL008426	(R)‐2,6‐dihydroxy‐2‐(4‐hydroxybenzyl)‐4‐methoxybenzofuran‐3(2H)‐one	49.8104	0.25786	Acacia catechu
5	MOL008432	Fisetinidol	49.63751	0.21042	Acacia catechu
6	MOL008421	Cis‐dihydro quercetin	47.73094	0.26823	Acacia catechu
7	MOL000098	Quercetin	46.43335	0.27525	Acacia catechu
8	MOL008430	5,7‐dihydroxy‐2‐(4‐hydroxyphenyl)‐3‐methyl‐chromone	45.04639	0.23585	Acacia catechu
9	MOL000422	Kaempferol	41.88225	0.24066	Acacia catechu
10	MOL008471	Isorhyncophylline	47.31	0.57	Acacia catechu
11	MOL002914	Eriodyctiol (flavanone)	41.35043	0.2436	Acacia catechu/ Scutellariae radix
12	MOL000073	Ent‐Epicatechin	48.95984	0.24162	Acacia catechu/ Scutellariae radix
13	MOL002934	NEOBAICALEIN	104.3446	0.43917	Scutellariae radix
14	MOL002932	Panicolin	76.25705	0.2915	Scutellariae radix
15	MOL012246	5,7,4′‐trihydroxy‐8‐methoxyflavanone	74.23522	0.26479	Scutellariae radix
16	MOL002927	Skullcapflavone II	69.51043	0.4379	Scutellariae radix
17	MOL002937	DIHYDROOROXYLIN	66.06174	0.23057	Scutellariae radix
18	MOL000228	(2R)‐7‐hydroxy‐5‐methoxy‐2‐phenylchroman‐4‐one	55.23317	0.20163	Scutellariae radix
19	MOL002915	Salvigenin	49.06593	0.33279	Scutellariae radix
20	MOL002917	5,2′,6′‐Trihydroxy‐7,8‐dimethoxyflavone	45.04743	0.33057	Scutellariae radix
21	MOL008206	Moslosooflavone	44.08796	0.25331	Scutellariae radix
22	MOL000449	Stigmasterol	43.82985	0.75665	Scutellariae radix
23	MOL001490	Bis[(2S)‐2‐ethylhexyl] benzene‐1,2‐dicarboxylate	43.59333	0.34531	Scutellariae radix
24	MOL002879	Diop	43.59333	0.39247	Scutellariae radix
25	MOL002897	Epiberberine	43.09233	0.7761	Scutellariae radix
26	MOL002928	Oroxylin a	41.36757	0.23233	Scutellariae radix
27	MOL002910	Carthamidin	41.15096	0.24189	Scutellariae radix
28	MOL002913	Dihydrobaicalin_qt	40.03778	0.20722	Scutellariae radix
29	MOL000525	Norwogonin	39.40397	0.20723	Scutellariae radix
30	MOL010415	11,13‐Eicosadienoic acid, methyl ester	39.27534	0.2289	Scutellariae radix
31	MOL012266	Rivularin	37.94023	0.3663	Scutellariae radix
32	MOL002925	5,7,2′,6′‐Tetrahydroxyflavone	37.01349	0.24382	Scutellariae radix
33	MOL000358	Beta‐sitosterol	36.91391	0.75123	Scutellariae radix
34	MOL000359	Sitosterol	36.91391	0.7512	Scutellariae radix
35	MOL012245	5,7,4′‐trihydroxy‐6‐methoxyflavanone	36.62689	0.26833	Scutellariae radix
36	MOL002933	5,7,4′‐Trihydroxy‐8‐methoxyflavone	36.562	0.26666	Scutellariae radix
37	MOL001689	Acacetin	34.97357	0.24082	Scutellariae radix
38	MOL002909	5,7,2,5‐tetrahydroxy‐8,6‐dimethoxyflavone	33.81583	0.44739	Scutellariae radix
39	MOL002714	Baicalein	33.51892	0.20888	Scutellariae radix
40	MOL000552	5,2′‐Dihydroxy‐6,7,8‐trimethoxyflavone	31.71246	0.35462	Scutellariae radix
41	MOL000173	Wogonin	30.68457	0.22942	Scutellariae radix
42	MOL001458	Coptisine	30.67185	0.85647	Scutellariae radix

Abbreviations: DL, drug‐likeness; OB, oral bioavailability.

Collection of target information and construction of active component‐target network

From the TCMSP and SIB platform, 367 targets of the 42 active compounds were obtained, including 250 from AC and 117 from SR. After eliminating repeated targets, 288 targets and their abbreviations remained. The gene names corresponding to the proteins were found in UniProt. The active ingredient‐target network was constructed using the network analysis software Cytoscape (version 3.7.1), and results are depicted in Figure 2. The octagon node represented the active ingredients of AC‐SR formula, and the triangle node represented the target genes. Edges represent interaction between compounds and targets. The size of the shape represented the degree of node association, the more connected the edges are, and a higher degree value is obtained. We list the top 8 components in Table 2 according to degree value between the components and the targets. It included MOL000098‐quercetin, MOL54758660‐Fisetin(1‐), MOL000422‐kaempferol, MOL000173‐Wogonin, MOL000358‐Beta‐sitosterol, MOL002714‐Baicalein, MOL002927‐Skullcapflavone II, and MOL000449‐Stigmasterol, with 143, 100, 61, 45, 38, 36, 33, and 31°, respectively.

FIGURE 2

Herb‐compound‐target network of AC‐SR formula (The ellipse nodes are composed of all the herbs of AC‐SR formula, which are surrounded with their particular compounds. The octagon nodes represent the compounds of AC‐SR formula. The Triangle nodes, arranged into a rectangular matrix, represent the relative gene targets of AC‐SR formula)

TABLE 2

Basic information for the top eight scored compounds

Mol ID	Compound	Pubchem CID	Molecular formula	Degree	Structure
MOL000098	Quercetin	5,280,343	C15H10O7	143
/	Fisetin(1‐)	54,758,660	C15H9O6	100
MOL000422	Kaempferol	5,280,863	C15H10O6	61
MOL000173	Wogonin	5,281,703	C16H12O5	45
MOL000358	Beta‐sitosterol	222,284	C29H50O	38
MOL002714	Baicalein	5,281,605	C15H10O5	36
MOL002927	Skullcapflavone II	124,211	C19H18O8	33
MOL000449	Stigmasterol	5,280,794	C29H48O	31

Potential targets of AC‐SR formula in the treatment of COVID‐19

A search of the GeneCards, DisGeNET, TTD, and CTD databases identified 7,659 target genes linked with COVID‐19. Venn diagrams indicated a total of 209 AC‐intersection (Figure 3a) and 99 SR‐intersection (Figure 3b) targets against COVID‐19. Then, a PPI network was constructed to integrate the targets to obtain the intersection. The targets are represented by the circle nodes, while interaction between targets is represented by edges. The degree of node linkage is symbolized by the size and depth of the circle, the more linked the edges are, the greater the degree value. Finally, we obtained 7 key targets from AC and 6 key targets from SR. After eliminating repeated targets, TP53, JUN, ESR1, MAPK1, Akt1, HSP90AA1, TNF, IL‐6, SRC, and RELA were considered to be the hub genes.

FIGURE 3

Targets of AC and SR against COVID‐19. (a) Venn diagram and PPI network showed the 209 targets of AC against COVID‐19. (b) Venn diagram and PPI network exhibited the 99 targets of SR against COVID‐19

GO functional enrichment and KEGG pathway enrichment analyses

To further analyze the target genes, 209 AC‐intersection and 99 SR‐intersection targets against COVID‐19 were implemented using Org.Hs.eg.db and ClusterProfiler packages in R for GO analysis and KEGG pathway analysis, respectively. The results showed the top 8 significantly enriched terms in biological processes (BP), cellular components (CC), and molecular functions (MF) (Figure 4a, p < 0.05). BP mainly included response to drug, cellular response to chemical stress, response to oxidative stress, etc. CC mainly included membrane raft, membrane microdomain, protein kinase complex, etc. MF mainly included DNA‐binding transcription factor binding, carbonate dehydratase activity, RNA polymerase II‐specific DNA‐binding transcription factor binding, etc. Similarly, the GO function histogram of SR in COVID‐19 treatment was shown in Figure 4c.

FIGURE 4

Biological processes and molecular pathways associated with core targets of AC (a, b) and SR (c, d) against COVID‐19. Biological processes [from GO analysis] were presented the top 8 BP (biological processes), CC (cellular components) and MF (molecular functions) by bar diagrams with count algorithms and p‐adjust values. Molecular pathways (from KEGG analysis) were presented the top 30 ranking pathways by bar diagrams based on ‐log10 (p‐adjust values)

The first 20 enriched KEGG pathways of AC in the treatment of COVID‐19 were shown in Figure 4b (p < 0.05). The Y‐axis represents the name of the pathway, the X‐axis represents the ratio of targeted genes to background genes, the size of the dot represents the number of genes concentrated on the modified pathway, and the color of the dot represents the significance of enrichment. AC‐KEGG pathway analysis targets mainly involved lipid and atherosclerosis, human cytomegalovirus infection, fluid shear stress and atherosclerosis, IL‐17, TNF, etc. signaling pathways. And SR‐related pathways were associated with lipid and atherosclerosis, human cytomegalovirus infection, PI3K‐Akt, AGE‐RAGE, IL‐17, TNF, etc. signaling pathways (Figure 4d), indicating that AC and SR may have synergistic reaction through multiple targets and multiple pathways in the treatment of COVID‐19. Figure 5 showed the relevant targets in the key inflammatory pathways of AC‐SR formula.

FIGURE 5

KEGG pathway enrichment map. (a) TNF signaling pathway, (b) PI3K/Akt signaling pathway, (c) IL‐17 signaling pathway. Red rectangles represent key targets

Molecular docking

We selected 8 core target proteins (TP53, JUN, ESR1, MAPK1, Akt1, HSP90AA1, TNF, IL‐6) in PPI and 2 of most important targets for treating COVID‐19 (3CLpro, ACE2) as the protein receptors via the AutoDock Vina software. We also added (+)‐catechin and (−)‐Epicatechin, the main active ingredients extracted from AC (Khare, 2007; Xu et al., 2015), and three positive control compounds (ritonavir, Nirmatrelvir and SSAA09E2) as ligands for molecular docking verification. The affinity between these compounds and the targets was lower than −5.0 kcal/mol, indicating that the core active compounds had a good binding activity with the main target. As shown in Table 3, 43 pairs of compound‐target interactions were of good binding affinity. Beta‐sitosterol, baicalein, stigmasterol showed strong affinity for Akt1, ESR1, HSP90AA1 with a binding free energy of less than −7.0 kcal/mol, respectively. The 10 compounds of AC‐SR formula were docking with 3CLpro and ACE2 in Table 4. Compared with the positive control drugs for 3CLpro (ritonavir, Nirmatrelvir), Baicalein, (+)‐catechin, and Fisetin(1‐) showed superior affinity with a binding free energy of less than −7.4 kcal/mol. Meanwhile, except for Skullcapflavone II and Wogonin, all compounds exhibited more excellent binding affinity for ACE2 with a binding free energy of less than −8.6 kcal/mol, compared with SSAA09E2. The docking position of the 8 typical compound‐target interactions is shown in Figure 6. The 2D image of molecular docking depicts the binding mode between the core target protein and the compound as well as the interaction with the surrounding amino acid residues. The main forces between ligand and protein are hydrophobic force and hydrogen bond, which both are chemical bonds with strong binding force (Fan et al., 2021). For instance, Beta‐sitosterol bounded to Akt1, had hydrophobic interactions with Phe150, Tyr184, Pro187, Ile156, Met176, Phe230, Met166, Ala162, His152, Leu110, and Pro147 (Figure 6a). When binding to ACE2, Stigmasterol formed hydrogen bond interactions with residues Tyr202, formed hydrophobic interactions with His401 (Figure 6d). Baicalein formed three hydrogen bonds with the amino acid residues and had hydrophobic interactions with Met49 and Gly143 as well, which make it form a stable complex with 3CLpro (Figure 6h).

TABLE 3

The affinity of the putative compounds with 8 core targets

Herb	Mol ID	Compound	Affinity value (kcal/Mol)
			ESR1 (1gwq)	JUN (1a02)	TNF (6 m95)	HSP90AA1(3o0i)	TP53 (4agq)	MAPK1 (4fv4)	AKT1 (4gah)	IL‐6 (1il6)
AC	MOL000098	Quercetin	−6.2	−4.4	−5.1	−6.5	−5.1	−4. 9	−5.1	−5.5
AC	/	Fisetin(1‐)	−6.9	−5.5	−5.4	−6.5	−4.7	−4.4	−5. 5	−5.4
AC	MOL000422	Kaempferol	−6.3	−6.4	−5.5	−6.2	−5.1	−5.4	−4.6	−5.1
SR	MOL000173	Wogonin	−2.9	−5.9	−3.0	−4.9	−2.5	−3.1	−2.5	−4.8
SR	MOL000358	Beta‐sitosterol	−6.1	−5.7	−5.1	−4.9	−5.1	−5.1	−7.5	−5.4
SR	MOL002714	Baicalein	−7.3	−4.9	−6.0	−6.7	−6.6	−5.2	−6.8	−5.2
SR	MOL002927	Skullcapflavone II	−5.0	−4.1	−4.3	−5.4	−4.6	−4.7	−4.9	−4.1
SR	MOL000449	Stigmasterol	−6.5	−5.9	−5.4	−7.2	−4.5	−5.4	−6.6	−5.6
/	Positive drug	Ritonavir	−5.2	−4.0	−5.0	−6.1	−4.9	−5.2	−5.4	−5.4
/	Positive drug	Nirmatrelvir	−5.1	−4.5	−5.8	−5.7	−5.1	−5.0	−5.3	−4.8

The binding energy ≤0 kcal/mol indicated that the compound could bind and interact with the target, whereas the binding energy <−5 kcal/mol indicated a very strong binding force.

TABLE 4

The binding energies of 10 compounds and positive drugs to 3CLpro and ACE2

Pubchem CID	Molecule name	Molecular formula	MW (g/Mol)	SARS‐CoV‐2 3CLpro (6lu7) docking score (kcal/Mol)	ACE2 (1r4l) docking score (kcal/Mol)
5,280,343	Quercetin	C15H10O7	302.25	−7.2	−9.1
54,758,660	Fisetin(1‐)	C15H9O6	285.23	−7.5	−8.8
5,280,863	Kaempferol	C15H10O6	286.25	−7.3	−8.8
5,281,703	Wogonin	C16H12O5	284.26	−6.7	−8.3
222,284	Beta‐sitosterol	C29H50O	414.79	−6.8	−9.5
5,281,605	Baicalein	C15H10O5	270.25	−7.8	−9.1
124,211	Skullcapflavone II	C19H18O8	374.30	−7.0	−8.2
5,280,794	Stigmasterol	C29H48O	412.77	−6.8	−9.8
9,064	(+)‐Catechin	C15H14O6	290.27	−7.8	−9.0
72,276	(−)‐Epicatechin	C15H14O6	290.27	−7.1	−8.8
392,622	Ritonavir (positive drug)	C37H48N6O5S2	720.96	−6.7	/
155,903,259	Nirmatrelvir (positive drug)	C23H32F3N5O4	499.50	−7.4	/
2,738,575	SSAA09E2 (positive drug)	C16H20N4O2	300.36	/	−8.6

The binding energy ≤0 kcal/mol indicated that the compound could bind and interact with the target, whereas the binding energy <−5 kcal/mol indicated a very strong binding force.

FIGURE 6

Molecular models of the selected compounds binding to the target proteins. (a) The docking mode and interactions between Beta‐sitosterol and AKT1(4gah), (b) Baicalein‐ESR1(1gwq), (c) Stigmasterol‐HSP90AA1(3o0i), (d) Stigmasterol‐ACE2(1r42), (e) Kaempferol‐ACE2(1r42), (f) Quercetin‐ACE2(1r42), (g) Stigmasterol‐3CLpro(6lu7), (h) Baicalein‐3CLpro(6lu7)

Molecular dynamic simulation

Based on the interactions with the binding pockets and binding energy calculations (Figure 6, Table 4), ACE2‐Quercetin and 3CL‐Baicalein complexes were selected to run MD simulations. The complex stability was predicted based on the RMSD and RMSF calculation during the entire run of MD trajectories. Figure 7a displayed the RMSD value of these two complexes over 200 ns run, and each trajectory was found with an average RMSD of 2.21–2.99 Å for ACE2‐Quercetin complex, 1.25–2.17 Å for 3CL‐Baicalein complex, respectively. The complex was stable over the entire 200 ns run. For gaining more insights regarding the stability of the complex binding site, the per‐residue RMSF was estimated for each ligand‐bound protein. The RMSF value was suggested about the very low fluctuations of molecule form the protein over entire MD trajectories (Figure 7b, c). Further, heatmap of hydrogen bonding interactions for ACE2‐Quercetin and 3CL‐Baicalein demonstrated the hydrogen‐bonding pattern (red blocks) observed during 200 ns simulation in both 2 protein‐ligand complexes (Figure 7d, e). These findings all showed that a stable conformation has been achieved in the process of MD simulation.

FIGURE 7

Molecular dynamics simulations. (a) The RMSD plot of ACE2‐Quercetin and 3CLpro‐Baicalein; The RMSF of ACE2‐Quercetin (b) and 3CLpro‐Baicalein (c); Heatmap of hydrogen bonding interactions for ACE2‐Quercetin (d) and 3CLpro‐Baicalein (e), red indicates hydrogen‐bonded, whereas green indicates non‐hydrogen‐bonded

FRET‐based assay for the SARS‐CoV‐2 3CLpro and ACE2 enzyme activity inhibition

The in silico docking study and molecular dynamic simulation assay both indicated that the active compounds have certain affinity with 3CLpro and ACE2, the key proteases during SARS‐CoV‐2 replication. The anti‐COVID‐19 property of these compounds may possibly due to its ability to inhibit the enzymatic activity of 3CLpro and ACE2. In this regard, the in vitro enzymatic inhibitory assays based on fluorescence resonance energy transfer (FRET) were applied to determine the median inhibitory concentration (IC50) values. The IC50 for the positive control Ebselen and MLN‐4760 have been calculated to be 0.67 μM and 7.5 nM, respectively. Among the top three scored compounds in Table 4, baicalein has the best inhibitory effect on 3CLpro activity with IC50 of 11.3 μM, followed by fisetin(1‐) and (+)‐catechin with IC50 of 23.8 and 44.1 μM, respectively (Figure 8). A concentration‐dependent inhibition of baicalein, quercetin and (+)‐catechin against ACE2 was also observed with the IC50 of 138.2, 141.3 and 348.4 μM, respectively. In addition, IC50 values for both beta‐sitosterol and stigmasterol were determined to be > 1,000 μM.

FIGURE 8

Dose–response inhibition of SARS‐CoV‐2 3CLpro and ACE2 activity by the selected compounds. For each compound, at least three independent experiments were performed for the determination (n = 3). Non‐linear regression (curve fit) with log (inhibitior) vs. response‐Variable slope was used to calculate the IC50 values

DISCUSSION

The COVID‐19 outbreak has sparked worldwide alarm as an emerging infectious disease. The practice of China in controlling this outbreak has demonstrated the clinical responses and superiorities of TCM (Jiang et al., 2021). Among the drugs recommended by the government and doctors, SR is one of the most frequently utilized herbs (60.84%) for treating COVID‐19 (Luo et al., 2020). Indeed, the Chinese medicinal formula AC‐SR, also famous as Huang‐qin Er‐cha Decoction, has been widely used for treating cough, phlegm and fever which caused by pulmonary infection in clinical practice for several centuries. It was reported that the pathophysiology of COVID‐19 has been linked to cytokine storms and consequent immunogenic damage. The elevated level of proinflammatory cytokines, such as IL‐6, IL‐1, and IFN‐γ, are related to the influx of leukocytes, which further accelerates the local inflammatory response in the lungs (Zhang et al., 2019). At this point, the disease may quickly progress to severe sickness, manifesting as ARDS, acute lung injury (ALI), multiple organ failure, and septic shock (Kim & Hong, 2016; Yang et al., 2020). Interestingly, our previous studies have demonstrated that AC‐SR formula exhibits a strong anti‐inflammatory effect against LPS‐induced ALI. Thus, to reveal the integrative pharmacological mechanism of AC‐SR formula against COVID‐19, we applied a unique systems pharmacology strategy which combining network pharmacology molecular, docking analysis, and MD simulations in the present study.

Eight key candidate components (quercetin, fisetin(1‐), kaempferol, wogonin, beta‐sitosterol, baicalein, skullcapflavone II, and stigmasterol) from AC‐SR formula were screened out based on the degree value of correlation between components and targets. The results of PPI network revealed that TP53, JUN, ESR1, MAPK1, Akt1, HSP90AA1, TNF, IL‐6 were considered to be hub genes. These proteins are mainly involved and play important roles in inflammation and immune regulation. JUN, for instance, is an immediate‐early gene that plays a crucial role in inflammatory responses (Fahmy et al., 2006). Akt1 activation promotes cell proliferation while inhibiting cell apoptosis. It is a key participant in the immune inflammatory mechanism of COVID‐19 (Xia et al., 2020). Furthermore, IL‐6 and TNF, which are critical components of the body’s immunomodulatory and proinflammatory effects, have also been well documented in the literature (Yi et al., 2020).

Notably enriched GO biological processes included response to oxidative stress and response to molecule of bacterial origin. AC‐KEGG enrichment analysis showed that the key targets were mainly concentrated in lipid and atherosclerosis, human cytomegalovirus infection, fluid shear stress, and atherosclerosis, etc. SR‐KEGG enrichment analysis was mainly concentrated in lipid and atherosclerosis, human cytomegalovirus infection, hepatitis B, etc. Moreover, according to our KEGG pathway analysis, the pharmacological mechanisms of AC‐SR formula against COVID‐19 involved specific modulations of immunological responses, such as human immunodeficiency virus 1 infection, human T‐cell leukemia virus 1 infection, and Human cytomegalovirus infection.

The deteriorated clinical presentation of patients with COVID‐19 is mainly associated with cytokine release syndrome (Zhang et al., 2020). It has been shown that the consequent inflammatory waterfall factors of host infection are critical for defense against and treatment of new coronavirus infections (Rokni, Ghasemi, & Tavakoli, 2020). Our data implicated TNF signaling pathway, PI3K/Akt signaling pathway, and IL‐17 signaling as the key inflammatory signaling pathways of AC‐SR formula anti‐inflammatory function to treat COVID‐19. TNF is a potent cytokine exerting critical functions in the activation and regulation of immune and inflammatory responses (Shivappa et al., 2017). TNF and TNF receptor 2 signaling elicited leukocyte recruitment, activation, and survival of host cells after coronavirus infection (Cheng et al., 2021). The PI3K/Akt signaling pathway regulated the activation of inflammatory response cells and the release of inflammatory transmitters. Targeting PI3K and Akt during the early phases of the immune response may improve effector function and suppress suppressor function, thereby assisting in the elimination of the infection before it progresses to immunological dysregulation. When patients develop uncontrolled immune responses, targeting mTOR can also be used to suppress the cytokine storm, inhibit neutrophil recruitment, and enhance suppressor regulatory T cells (Abu‐Eid & Ward, 2021). IL‐17 family proteins not only mediated the cytokine storm after SARS‐CoV‐2 infection (Cafarotti, 2020) but also regulated the innate immune responses (Ryzhakov et al., 2011). IL‐17 signaling played an underlying immunopathological role to manage COVID‐19 patients, particularly those presenting with cytokine storm syndrome (Shibabaw, 2020).

ACE2 and 3CLpro have been considered as promising COVID‐19 drug targets, which play a key role in the replication cycle of the virus. Researches showed that ACE2 and 3CLpro on host epithelial cells affected by its S‐protein are considered to be the core targets for inhibiting coronavirus proliferation (Hall & Ji, 2020; Menachery et al., 2015; Wu et al., 2020). In molecular docking study, the results showed Beta‐sitosterol, baicalein, stigmasterol showed strong affinity for Akt1, ESR1, HSP90AA1, respectively. Compared with the positive control drugs for 3CLpro and ACE2, Baicalein, (+)‐catechin, Fisetin(1‐), beta‐sitosterol, stigmasterol, and quercetin showed superior or similar affinity with a lower binding free energy. Beta‐sitosterol and stigmasterol are the most frequently used chemical components possibly related to the antiviral signing pathway (Luo et al., 2020). Unfortunately, IC50 values for both beta‐sitosterol and stigmasterol were determined to be greater than 1,000 μM in the ACE2 activity inhibition assay, probably due to their poor solubility in water and alcohol, thus further evaluation of the two compounds as potential SARS‐CoV‐2 ACE2 inhibitors may be warranted. Particularly, we identified that baicalein had good affinity for ACE2 and 3CLpro, mainly through hydrogen bonds and hydrophobic interactions, with measured IC50 values against SARS‐CoV‐2 3CLpro and ACE2 to be 11.3 and 138.2 μM, respectively, which indicate that baicalein may directly inhibit SARS‐CoV‐2 replication and transcription. Du et al. (2022) also confirmed that baicalein has good binding activity with ACE2 and that hydrogen bonding plays a key role in the recognition and stability of the active ingredients and proteins. In addition, baicalein is the first identified non‐covalent, non‐peptidomimetic inhibitors of SARS‐CoV‐2 3CLpro (Liu et al., 2021; Su et al., 2020). Notably, We found for the first time that fisetin(1‐) might have the potential against the novel coronavirus via binding to SARS‐CoV‐2 3CLpro and ACE2 in silico. This was further confirmed by the enzyme activity inhibition assays. Fisetin(1‐) revealed the prominent inhibitory activity against SARS‐CoV‐2 3CLpro with the IC50 of 23.8 μM. We also found that quercetin and kaempferol, which have multiple pharmacological activities such as anti‐inflammatory, immunomodulatory and antiviral (Carullo et al., 2016; Chiow, Phoon, Putti, Tan, & Chew, 2016; Huang, Bai, He, Xie, & Zhou, 2020; Khazdair, Anaeigoudari, & Agbor, 2021; Mlcek, Jurikova, Skrovankova, & Sochor, 2016), plays crucial roles against COVID‐19 via binding to ESR1, HSP90AA1, TP53, TNF, ACE2, and 3CLpro. These findings were in line with several previous studies (Huang et al., 2021; Qiu et al., 2020; Xia et al., 2020). Catechin, the main antioxidant, anti‐inflammatory, and chemoprotective ingredient of AC, exhibits excellent affinity for ACE2 and 3CLpro in the docking test as well. Several previous work have reported that catechin exhibit strong in silico activity against the wild strain of SARS‐CoV‐2 (Mhatre, Gurav, Shah, & Patravale, 2021), presumably leading to therapeutic efficacies via inhibition of ACE2 and could be potentially explored as a multitargeted agent against COVID‐19 (Jena, Kanungo, Nayak, Chainy, & Dandapat, 2021; Mishra et al., 2021). Considering this, quercetin, fisetin(1‐), kaempferol, catechin, beta‐sitosterol, baicalein, and stigmasterol could be the potential candidates of AC‐SR formula against COVID‐19. In addition, in order to validate the results of molecular docking as well as evaluate the stability and interaction of ACE2 and 3CLpro with ligands, ACE2‐quercetin and 3CLpro‐baicalein complexes were chosen for molecular dynamic simulation test. The binding sites of small molecules did not change considerably after a 200 ns dynamics simulation, suggesting that the results of molecular docking were reliable and the binding of small molecules remained relatively stable. The two complexes also exhibited stable RMSD and RMSF in the simulated test, which demonstrated that quercetin and baicalein may effectively activate the biological pathway in ACE2 and 3CLpro without affecting the conformation of the active site. The amount of hydrogen bonds as well as hydrophobic bond, on the other hand, influences the stability of protein‐ligand complexes. The aggregates formed by 3CLpro‐baicalein and ACE2‐quercetin in the MD test contained 3–5 hydrogen bonds and hydrophobic bond, which suggested that the strong interactions of the key residues (most especially Tyr196, Ser144, Cys145, Leu141, Lys562 and Gln96) with quercetin and baicalein would significantly impede SARS‐CoV‐2 3CLpro and ACE2 dimerization and substrate binding.

CONCLUSIONS AND PERSPECTIVES

In summary, the present study dissected the possible molecular targets and mechanisms of Chinese medicinal formula AC‐‐SR against COVID‐19 from the point of view of a systems pharmacology strategy through network pharmacology approach, molecular docking and molecular dynamic simulation methods. After screening and analysis, 8 active ingredients and 10 core target proteins in AC‐SR formula were obtained, involving multiple pathways. AC‐SR formula might suppress COVID‐19 through their combined antioxidative, antiviral and anti‐inflammatory effects, along with immune system activation. The present work provides a rapid and accurate method for facilitate the screening and dissecting of the complex system of TCM and natural medicine. Based on the identified functional processes and pharmacological mechanisms, the Chinese medicinal formulas as well as their ingredients, could be applied to the development of anti‐SARS‐CoV‐2 medications in a more effective and targeted manner. However, for the identification and mechanism study of novel COVID‐19 medication combination, we believe that the database‐based systems pharmacology should be supplemented by experimental evidence, which will not only decrease the duration of drug development in crisis, but will also assure its dependability. Therefore, in light of our findings, related and more in‐depth studies, such as direct evidence for the mechanism of AC‐SR formula in a SARS‐CoV‐2 infection experiment model, remains urgently warranted. The further verification of our present study will be on the way.

ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China (No. 81703737), Key Research & Development Project of Shaanxi Province (Nos. 2020SF‐209, 2021SF‐152 and 2022SF‐253), Basic Research plan of Natural Science of Shannxi Province (No. 2021JQ‐887), and Key Science and Technology Program for Social Development of Shaanxi Province (No. 2016SF‐395).

CONFLICT OF INTEREST

The authors declare no conflict of interest.