A Network Pharmacological Approach to Explore the Mechanisms of TongXieYaoFang in Inflammatory Bowel Disease

Inflammatory bowel disease (IBD) is a chronic inflammatory disorder of the intestine, demonstrating an increasing incidence every year. TongXieYaoFang (TXYF) has been used widely in China as a complementary therapy to relieve the symptoms of IBD for hundreds of years. In the present research, a network pharmacology-based approach was used to systematically explore the intrinsic mechanisms of TXYF in IBD at the molecular level. Network pharmacology-based methods, which mainly included database mining, screening of bioactive compounds, target prediction, collection of IBD-related targets, gene enrichment analysis, network construction, and molecular docking, were employed in the present study. Network analysis revealed a total of 108 potential targets derived from 22 component compounds of TXYF, among which 34 targets were common with the IBD-related targets. In the protein-protein interaction (PPI) network, 10 key targets were identified. The gene enrichment analysis suggested that anti-inflammatory processes, such as NF-kappa B signaling pathway and Toll-like receptor signaling pathway, could be the core processes involved in the action of TXYF in IBD. Molecular docking results revealed that three compounds present in TXYF exhibited strong binding affinity for PTGS2. The present study provides novel insights into the molecular mechanisms and network approaches of TXYF action in IBD from a systemic perspective. The potential targets and pathways identified in the present study would assist in further research on the clinical application of TXYF in IBD therapy.


INTRODUCTION
Inflammatory bowel disease (IBD) is a chronic diffuse inflammatory disorder of the intestine, which is characterized by abdominal pain and diarrhea in the intestinal mucosa (Sairenji, Collins, Evans, 2017;Zhang, Li, 2014). In the last few decades, the incidence of IBD has increased exponentially worldwide, with a higher prevalence in North America and Europe compared to Asia and Africa (Sairenji, Collins, Evans, 2017). Patients with IBD present accelerated inflammatory responses and a significantly increased risk of colorectal cancer (Axelrad, Lichtiger, Yajnik, 2016). Although immunosuppressive and anti-inflammatory drugs are recommended for IBD patients to prevent or delay the expected inflammatory events, the majority of the patients do not exhibit any significant benefits in terms of full recovery (Feuerstein et al., 2018;Mitrev et al., 2017;Peyrin-Biroulet et al., 2010). Therefore, there is an urgent requirement for a feasible and effective therapeutic strategy capable of alleviating inflammation for the treatment of IBD. Traditional Chinese medicine (TCM), which has been traditionally used widely in the treatment of inflammatory diseases due to its multiple pharmacological effects, has emerged as a safe and effective complementary therapeutic strategy for the treatment of IBD patients. TongXieYaoFang (TXYF), also known as BaiZhuShaoYao San, was first documented in "DanXi-Xin-Fa" (1281 A.D.). The TCM formula of TXYF comprises four medicinal herbs, namely, Atractylodes Macrocephala Koidz (Baizhu), Paeoniae Radix Alba (Shaoyao), Citrus Reticulata (Chenpi), and Saposhnikoviae Radix (Fangfeng). It is reported that TXYF treatment is effective against abdominal pain and diarrhea (Xu et al., 2017). Moreover, TCM theory as well as practice has demonstrated that prescription compatibility offers the advantages of fewer side effects and better efficacy over the use of a single herb (He et al., 2015;Zhao et al., 2015). Recent pharmacological studies have indicated that TXYF has a significantly higher clinical effective rate and decreases the adverse events related to Diarrhea Predominant Irritable Bowel Syndrome (Bian et al., 2006;Zhou, Han, He, 2019). Furthermore, Cai et al. (2019) observed that TXYF exerted an inhibitory effect on the inflammatory response elicited in UC rats, as indicated by the results of UHPLC/Q-TOF-MS/MS analysis. However, the functional mechanism underlying the role of TXYF in the inflammatory response during IBD remains unknown largely.
Scholars working on deciphering the potential molecular mechanisms of TCM encounter great challenges due to the multi-component and multi-target mode of action of TCM. Network-based drug discovery is a fairly recent novel strategy for analyzing complex and multilayered networks using biological, pharmacological, and pharmacodynamics data and information (Hopkins, 2008). In the present study, network pharmacology analysis was used to explore the potential pharmacological mechanisms of action of TXYF in relieving IBD ( Figure 1). The results revealed a potential mechanism that could underlie the therapeutic effects of TXYF in IBD, thereby contributing to the development of novel drugs targeting IBD.

Active compounds of TXYF
The chemical components of the four herbs of TXYF were obtained from TCMSP (http://tcmspw. com/tcmsp.php) (Ru et al., 2014), a systematic pharmacological platform designed for comprehensive drug screening and evaluation. The main pointers for evaluating ADME (absorption, distribution, metabolism, and excretion) were oral bioavailability (OB) and drug-likeness (DL). The screening criteria weres OB ≥ 30% and DL ≥ 0.18. Among the 376 compounds of TXYF, 46 compounds were identified based on the criteria (Table I). Next, the standard names and 2-dimensional structures of these compounds were obtained using the PubChem CIDs (https://pubchem. ncbi.nlm.nih.gov/) (Kim et al., 2016). Finally, due to the lack of precise structural information for all compounds, only 22 compounds were retained for subsequent analysis. Informed consent was obtained from each patient prior to commencing the experiments.

Screening of target genes
DrugBank (https://go.drugbank.com/) (Probst, Reymond, 2018) and TCMSP databases were searched for retrieving the gene targets of the 22 active compounds screened-out in the previous step by selecting Homo sapiens as the species and using the screening threshold of ≥0.7. UniProt online tool (https://www.uniprot.org/) (UniProt Consortium, 2018) was employed to standardize the gene names of the screened-out targets and remove the duplicate data.

IBD target genes
Genes associated with IBD were obtained by searching the DisGeNET database (https://www. disgenet.org/) (Piñero et al., 2017) using the keyword "IBD", which revealed a total of 910 genes with a genedisease score of ≥0.1. In order to identify potential targets of TXYF for IBD treatment, the predicted targets of TXYF were compared with the obtained IBD-related genes, revealing 34 common genes as overlap targets.

Protein-protein interaction (PPI) network
A PPI network was constructed using the String database (https://string-db.org/cgi/input.pl) (Szklarczyk et al., 2017), which is used for forecasting proteinprotein interactio verlap targets were submitted to this database, with the species set to "Homo sapiens" and the threshold criterion set to confidence score > 0.4. The top ten proteins were selected to construct the core PPI subnetwork, mainly using the Cytohubba plug-in (Chin et al., 2014), based on degree that could reflect the importance of the proteins.

Network construction
The following networks were constructed: (1) a network of the positive compounds of TXYF and the target genes; (2) a network of the overlap targets between the active components of TXYF and the IBD-related targets; (3) the PPI network of the key targets. The network constructions were visualized using Cytoscape version 3.7.2 (https://cytoscape. org/) (Shannon et al., 2003). Nodes in the networks represented targets and compounds, while the edges indicated interactions.

Biological function analysis
In order to investigate the functional annotation of the targets and the pathways involved, the overlap targets of TXYF and IBD were subjected to Gene Ontology (GO) (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa, Goto, 2000) enrichment analyses, followed by calculation and evaluation in the DAVID database version 6.8 (https:// david.ncifcrf.gov/) (Huang da, Sherman, Lempicki, 2009). GO analysis results were visualized using the clusterProfiler package in R software (v3.5.3) (Yu et al., 2012). KEGG analysis was conducted using the ClueGO plug-in (Bindea et al., 2009) in Cytoscape. P < 0.05 was considered statistically significant.

Molecular docking
The molecular docking simulation was performed to value the interactions. Crystal structures of PTGS2 were retrieved from the RCSB Protein Data Bank (PDB) (https://www.rcsb.org/) (Goodsell et al., 2020), which is a database providing archive information of the 3D structures of proteins and complex assemblies. The structure of the Crystal Complex of Cyclooxygenase-2 with indomethacin-butyldiamine-dansyl conjugate (PDB ID: 6BL3) was selected for this analysis. Maestro 12.0 (Schrödinger, LLC, New York, NY) (Beard et al., 2013;Salam et al., 2014;Zhu et al., 2014) was employed for the molecular preparation and docking procedures. Briefly, the crystal structure selected from PDB was processed utilizing the ProPrep module in Maestro to delete the water molecules, assign bond orders, add hydrogen atoms, assign zero-order bonds to metals and disulfide bonds, and generate heat states using Epik. In addition, the 3D structures of the selected compounds (wogonin, 5-O-Methyivisamminol, and beta-sitosterol) were obtained from PubChem and subsequently prepared for docking by correcting the protonation using the LigPrep module in Maestro with default options. Further, the glide module in Maestro was employed to dock the compounds using default parameters. The location of the original protein-ligand was defined as the active site, where the small molecule ligand binding sites were present.

Compound target network analysis
In this analysis, 22 bioactive components of TXYF (BaiZhu, ShaoYao, ChenPi, and FangFeng) were identified as candidate compounds for subsequent analyses using TCMSP and PubChem databases. These compounds, with OB and DL filters and favorable pharmacokinetic profiles, were retained as the active ingredients of TXYF. The 2D structures of these 22 compounds are depicted in Figure 2, and the detailed information is presented in Table I. In order to identify the targets of these TXYF constituents, a compound/ target network was constructed using Cytoscape. As depicted in Figure 3, the constructed compound/ target network has 134 nodes (4 herbal nodes, 22 active ingredient nodes, and 108 target nodes), implying the multi-target characteristics of the identified compounds. Note that beta-sitosterol is the main compound in three (ShaoYao, ChenPi, and FangFeng) of the four herbs of TXYF. In addition, over 20 genes were the main molecular targets of five compounds (wogonin, betasitosterol, 5-O-Methyivisamminol, kaempferol, and nobiletin). Finally, nine genes (PTGS2, TGS1, SVN5A, ESR1, DPP4, NCOA2, PIK3CG, AR, and GABRA1), regulated by over 8 components, were selected as the key targets of TXYF.

Candidate targets of TXYF against IBD
A total of 374 potential targets of the 22 active constituents of TXYF were obtained from DrugBank and TCMSP. Among these, 108 candidate targets were selected for subsequent analysis after standardized processing using UniProt. Meanwhile, a total of 910 gene targets associated with IBD were retrieved from the DisGeNET database. An intersection of the target genes of TXYF with the IBD-related genes revealed 34 common targets ( Figure   4A-4B). The intersection network suggested that these 34 common/shared targets could be regulated by multiple compounds (19/22). Among the 19 compounds depicted in Figure 4C, eight compounds, namely, wogonin, betasitosterol, 5-O-Methyivisamminol, kaempferol, nobiletin, divaricatol, ledebouriellol, and naringenin, could regulate a minimum of six IBD-related genes. In addition, a total of nine targets, namely, PTGS2, PTGS1, ESR2, BCL2, ACHE, NOS2, ESR1, DPP4, and PPARG, could be regulated by a minimum of five compounds.

GO and KEGG pathway enrichment analyses
DAVID bioinformatics online application was used for performing the GO enrichment and KEGG pathway analysis of the identified 34 targets using Homo sapiens as the gene symbol. The results were visualized using the R package. The GO analysis revealed that the mechanism of TXYF in the treatment of IBD was mainly related to the following biological processes (BP), molecular functions (MF), and cellular components (CC): enzyme binding, gene transcription, oxidation-reduction, and cytokine activity within the nucleus, cytosol, and extracellular space ( Figure 6A-6C). The KEGG pathway enrichment analysis revealed that IBD-related targets and TXYF-targeted genes were involved in the AGE-RAGE signaling pathway, TNF signaling pathway, Toll-like receptor signaling pathway, IL-17 signaling pathway, and NF-kappa B signaling pathway ( Figure 6D). In order to further investigate the mechanism of TXYF in the treatment of IBD, a compound-targetpathway network was constructed using the 19 compounds, 34 targets, and the corresponding pathways (Figure 7). Among the 34 targets, a total of 15 genes, namely, TNF, CCL2, MAPK8, MMP1, PTGS2, MAPK14, PLA2G4A, MAPK1, ALOX5, CCND1, NOS3, IL6, VCAM1, TGFB1, and BCL2, were identified as bridge targets between TXYF compounds and pathways.

PPI network analysis
The PPI network represents large-scale biological processes, which are fundamental to understanding cellular organizations (Yu, Zheng, 2019). A total of 34 common targets, which were associated with both anti-IBD activity and TXYF components, were imported into the STRING database for the construction of the PPI network. The network contained 34 nodes and 239 edges ( Figure 5A). Subsequently, the Cytohubba plug-in in the Cytoscape software was employed for a module analysis of the 34 targets in the PPI network. The clustering module revealed ten genes that had specific biological significance, namely, PTGS2, IL6, MAPK8, TP53, TNF, CCL2, MAPK1, MAPK14, PPARG, and HMOX1 ( Figure 5B). Wenli You,Mingjuan Li,Aiting Di,Xin Li1,,Hairui Gao,Cuixia Qiao,Bin Yu,Gang Zhao

Molecular docking simulation analysis
On the basis of network pharmacology results in the present study and previous reports, PTGS2 was identified as the target to be used in the docking analysis. Table II presents the docking scores and glide scores of the TXYF hit compounds (wogonin, betasitosterol, and 5-O-Methyivisamminol) against the active sites of the identified target PTGS2, obtained using the Glide module in Schrodinger software. According to the "lock and key principle" and the "induced fit" theory (Koshland, 1958), the ligand molecules are appropriately bound within the protein pocket based on the complementarity between the electrostatic distribution across the protein and the strong interaction forces. Glide score reflects the van der Waals forces and the electrostatic contributions (Halgren et al., 2004). In general, the smaller the Glide score, the stronger is the binding force. Among all the hit compounds, wogonin presented the lowest docking and glide scores against PTGS2, indicating a higher affinity between wogonin and PTGS2. The interaction diagram for wogonin in the active site of the PTGS2 protein ( Figure 8A-8B) illustrates the stabilization of its interaction through the formation of π-π stacking with the key residues TYR385 and TRP387. On the other hand, the interaction diagram ( Figure 8C-8D) of beta-sitosterol, which exhibited a high affinity for the PTGS2 protein, illustrates the formation of a hydrogen bond with PTGS2 at the critical site of GLU524 residue. Moreover, the two hydrogen bonds in 5-O-Methyivisamminol contributed to the stabilization of its interaction with the residues LYS83 and GLU524 of the PTGS2 protein ( Figure 8E-8F). The ligand interaction was also stabilized through the formation of a π-π stacking bond with the TYR115 residue of the PTGS2 protein.  -sitosterol -5.195 -5.195 5-O-Methyivisamminol -5.067 -5.067

DISCUSSION
Unsatisfactory efficacy of drugs against IBD affects the therapeutic outcomes in IBD treatment. There remains much to be investigated in respect to the pharmacological efficacy of herbal medicines when exploring the therapeutic strategies involving low-toxicity and highly-targeted drugs against IBD. TXYF is a commonly prescribed formula for the treatment of IBD in China. A few previous studies have reported the mechanism of action of certain components of TXYF in relieving IBD (Bose, Kim, 2013;Ho, Kuo, 2014). In the present study, the active components and targets of BaiZhu, ShaoYao, ChenPi, and FangFeng herbs (TXYF) were obtained from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database. The structures of these compounds were obtained from PubChem. Meanwhile, the target genes of IBD were screened out from the Dis-geNET database. Subsequently, the overlap targets between TXYF and IBD were identified and subjected to gene enrichment analysis and protein-protein interaction (PPI) analysis.
A total of 8 compounds among those evaluated were identified as the core compounds that acted on a minimum of 6 targets -wogonin, beta-sitosterol, 5-O-Methyivisamminol, kaempferol, nobiletin, divaricatol, ledebouriellol, and naringenin ( Figure 4C). Wogonin has been previously reported to ameliorate colitis in vitro as well as in vivo Sun et al., 2015). Beta-sitosterol is a plant sterol, which has been suggested to improve intestinal inflammation in a T cell-dependent manner (te Velde et al., 2015). Kaempferol, nobiletin, and naringenin have been demonstrated to inhibit multiple inflammatory pathways, in a response against IBD, in a number of in vitro studies (Al-Rejaie et al., 2013;Bian et al., 2019;He et al., 2018). Accordingly, it was inferred that these core compounds present in TXYF could play an important role in treating IBD. However, besides the aforementioned compounds that have already been validated, the remaining ones are yet to be validated in terms of their effectiveness.
Two vital targets, PTGS2 and BCL2, were identified as the key targets for TXYF in the treatment of IBD, as presented in Figure 4C, Figure 5B, and Figure 7. PTGS2 is the rate-limiting enzyme in the reaction process of the production of prostaglandins (Cox et al., 2005). Prostaglandins are thought to be essential for the process of wound healing in the gastrointestinal tract (Wallace, 2001). Therefore, PTGS2 expression might be a protective response aimed at curing. In addition, PTGS2 is reportedly linked to pain burden in IBD patients (Grossi et al., 2020). On the other hand, the anti-apoptotic protein BCL2 is reported to be highly expressed in IBD, which suggests the dysregulation of the apoptosis mechanisms in IBD (Dudzinska et al., 2019). AChE activity is reported to be significantly lower in IBD patients (Maharshak et al., 2013), although further details of the mechanism remain unknown so far.
Among the potential pathways revealed in the present study, NF-kappa B signaling pathway, Toll-like receptor signaling pathway, IL-17 signaling pathway, TNF signaling pathway, VEGF signaling pathway, Fc epsilon RI signaling pathway, and GnRH signaling pathway had a high degree value. A few of these pathways constitute significant mechanisms of IBD and have been proven to be influenced by TXYF, as discussed widely in the existing literature (Cai et al., 2019;Globig et al., 2014;McDonnell et al., 2011;Onizawa et al., 2009;Wang et al., 2003;Xu et al., 2017). However, IL-17 signaling pathway, VEGF signaling pathway, and Fc epsilon RI signaling pathway have not been investigated to date, which implies that the investigation of the mechanisms of action of TXYF remains incomplete currently.
Wogonin was originally isolated from Scutellaria Baicalensis Georgi. Wogonin exhibits a variety of biological activities, such as anti-inflammatory, antimicrobial, and anti-tumor activities. Beta-sitosterol is present in Shaoyao, Chenpi, and Fangfeng, and exerts anti-inflammatory effects by regulating cytokine levels. Molecular docking simulations have demonstrated that the interaction between wogonin and PTGS2 has higher agonistic activity, suggesting PTGS2 as a potential target for TXYF acting against UC.
However, one may only speculate rather than confirm these mechanisms as most of the results remain to be validated experimentally, which would be the future direction of our research work. In addition, further research to decipher the exact toxicological effects of TXYF is warranted.

CONCLUSION
In the present study, the potential mechanisms of TXYF acting against IBD were evaluated systematically using the pharmacology network analysis. Among the IBD-related genes, a total of 34 genes were identified as candidate targets of TXYF constituents evaluated in the present study. The prediction results also revealed eight core compounds associated with IBD and nine hub targets, namely, PTGS2, PTGS1, ESR2, BCL2, ACHE, NOS2, ESR1, DPP4, and PPARG, associated with both TXYF components and IBD treatment. Furthermore, the GO and KEGG analyses indicated that the effects of TXYF on IBD might be associated with the vital targets PTGS2 and BCL2 involved in the regulatory inflammatory signaling pathway. The present study provides a basis for further exploration of the molecular targets of TXYF by applying the network pharmacology approach to develop drug treatment against IBD.

DATA AVAILABILITY
The data used to support the findings of this study are available from the corresponding author upon request.