Utilizing network pharmacology and molecular docking to explore the underlying mechanism of Guizhi Fuling Wan in treating endometriosis

Chemical ingredients database building

To collect the chemical ingredients of the five herbs contained in GZFLW, the Traditional Chinese Medicine Systems Pharmacology Database (Ru et al., 2014) (TCMSP, https://tcmspw.com/tcmsp.php), the Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine (Liu et al., 2016) (BATMAN-TCM, http://bionet.ncpsb.org/batman-tcm/index.php) and Chinese Academy of Sciences Chemistry Database (Wang et al., 2015) (CASC, http://www.organchem.csdb.cn/scdb/default.htm) were used, which are both bioinformatics analysis tools for the main components of TCM. Four hundred ninety-four herbal ingredients were screened out.

Active ingredients screening

According to the characteristics of absorption, distribution, metabolism, and excretion (ADME) of drug, oral bioavailability (OB), and drug-likeness (DL) were used as screening indexes. OB is the percentage of oral drugs absorbed into the bloodstream, which is a frequently-used pharmacokinetic parameter. The OB was obtained using the OBioavail1.1 software, which covers 805 different drug and drug-like molecules (Liu et al., 2013). OB is a crucial indicator to judge whether the active ingredient can become a feasible therapeutic molecule (Xu et al., 2012). DL is a parameter to measure the ADME of drug molecules, which could help optimize pharmacokinetic and pharmaceutical properties. The DL threshold is 0.18 (Liu et al., 2013), which depends on the average of Drugbank (https://www.drugbank.ca/). DL is often used to select the active ingredients with “drug-like” properties in TCM composition (Tao et al., 2013). The OB and DL indices of all the related ingredients are presented in the TCMSP. In this process, those ingredients with OB ≧30% and DL ≧0.18 were chosen as the potential effect components for the next step.

GZFLW targets prediction

The active compounds’ targets in GZFLW were obtained from TCMSP, BATMAN-TCM, and CASC database, with the species limited as “Homo sapiens”.

Endometriosis therapeutic targets database building

Related targets of EMs were screened out by the GeneCards Database (Stelzer et al., 2016) (https://www.genecards.org/). The GeneCards database is often used to predict genetic information related to human diseases. We searched the database with the keyword “endometriosis” to obtain targets. Finally, the targets of GZFLW active ingredients were matched with the therapeutic targets of EMs. The overlapping targets were chosen as the potential targets of GZFLW in the treatment of EMs. These targets were then uploaded to STRING Database (Szklarczyk et al., 2017) (http://string-db.org/) to obtain the interactions of the screened targets with the confidence ≧0.7 and the result of protein-protein interaction (PPI).

Network topological feature set definitions

We chose two parameters to evaluate the topological characteristics: “Degree” reflects the number of other nodes interacting with this node; “Betweenness Centrality” (BC) is measured by the percentage of all shortest paths. The nodes with high betweenness can significantly impact the net by controlling the information transmitted between other nodes. The node parameters are positively correlated with their topology importance in the network.

Molecular docking of hub targets and active ingredients

Molecular docking is a computational tool to predict the binding ability and connection type of proteins and ligands. It can calculate and predict the conformation and direction of ligands at active protein sites. AutoDock 1.5.6 (http://autodock.scripps.edu/) was a molecular docking software that could be used to dock the hub targets and active ingredients based on network pharmacology. In this docking process, the 3D structure of thirteen hub targets were retrieved from RCSB Protein Data Bank (PDB) (http://www.rcsb.org/): AKT1 (PDB ID: 4EKL, 6S9W), TNF (PDB ID: 2E7A, 2ZJC, 2AZ5), TP53 (PDB ID: 4MZI, 6FF9), VEGFA (PDB ID: 4WPB, 3QTK, 4QAF), MAPK1 (PDB ID: 4G6N, 1WZY), MMP9 (PDB ID: 6ESM, 4WZV), JUN (PDB ID: 5T01, 1T2K), MAPK8 (PDB ID: 3VUK, 4L7F), INS (PDB ID: 4AJX, 4CY7), EGF (PDB ID: 1JL9,), IL6 (PDB ID: 4O9H, 1ALU), PTGS2 (PDB ID: 5F19), and FOS (PDB ID: 1FOS); the 3D shapes of active compounds were provided from ZINC Database (http://zinc.docking.org/) and PubChem Database (https://pubchem.ncbi.nlm.nih.gov/). Binding energy was used as docking score to evaluate the protein-ligand binding potential of molecular docking. Among them, those results with value ≤ -5 were selected and considered to have moderate binding potential and tight combination.

Enrichment analysis

We used the DAVID (https://david.ncifcrf.gov/) (Huang, Sherman & Lempicki, 2008) for GO enrichment analysis, including biological process (BP), molecular function (MF) and cellular component (CC). KEGG enrichment analysis was performed using the Metascape Database to obtain potential target-pathways (Zhou et al., 2019) (http://metascape.org/).

Network construction

Networks, including the hub target-compound and target-pathway, were then constructed. All networks were built using Cytoscape 3.7.2 (http://www.cytoscape.org/) (Shannon et al., 2003), an open-source software platform for visualization and data analysis of complex networks.

Primary endometrial stromal cells culture

HEM15a cells were cultured in DMEM containing 10% FBS and 1% penicillin/streptomycin at 37 ^∘C in the cell incubator with a humid atmosphere containing 5% CO 2. After that, the vimentin staining of cells was identified by immunofluorescence.

MTT colorimetric assay

The hEM15a cells were seeded in a 96-well plate and divided into the control group and the drug group. After incubation at 37 ^∘C for 48 h and reacted with GZFLW of different concentrations, hEM15a cells were reacted with 10 μL MTT solution. Finally, 150 μL dimethylsulfoxide (DMSO) was added after reacting for 4 h. The optical absorbance was detected at 568 nm by a plate reader.

Western blotting

The hEM15a cells were splitted in RIPA buffer (Beyotime) added with 10 μl PMSF (Aladdin) and phosphatase inhibitors (Beyotime) to extract total proteins. After that, equivalent amounts of proteins were resolved by poly acrylamide gel electrophoresis (PAGE) and transferred to polyvinylidene fluoride (PVDF) membranes (Millipore, Massachusetts, USA). Then, TBST with 5% skim milk was adopted for blocking. Afterward, the membranes were reacted with primary antibodies against GAPDH (1:1000, AB-P-R 001, Xianzhi, Hangzhou, China), BAX (1:1000, Ab32503, Abcam, Britain), P53 (1:1000, 10442-1-AP, Sanying, Wuhan, China), and caspase3 (1:2000, Ab184787, Abcam), at 4 ^∘C overnight. Subsequently, the matched secondary antibody was added to the membranes.

Apoptosis analysis by flow cytometry

The apoptosis of hEM15a cells was evaluated by Annexin V-FITC/PI apoptosis kit (KaiJi, Nanjing, China). The collected hEM15a cells, divided into normal group and GZFLW group, were resuspended in 500 ul binding buffer and added with Annexin V-PI solution. Then, the hEM15a cells were incubated for 5∼15 min without light. The percentage of apoptotic cells was then immediately detected on a flow cytometry (Beckmancoulter).

Active compounds of GZFLW

By retrieving from the TCMSP, BATMAN-TCM, and CASC database, there were 494 related ingredients of GZFLW in total, and there were 220 (44.5%) of CT, 34 (6.9%) of PC, 55 (11.1%) of CM, 119 (24.1%) of RPR, 66 (13.3%) of PK. With OB ≧30% and DL ≧ 0.18 as indexes, 89 active ingredients were screened out (Table 1), and the Herbs-Compounds network was constructed as Fig. 2.

Target prediction analysis

In this process, we collected 284 targets of 89 active ingredients, and there were 58 in CT, 48 in PC, 185 in CM, 104 in RPR, and 62 in PK. Via the keyword of “endometriosis,” the therapeutic targets of EMs were obtained from GeneCard Database with a total of 1350. One hundred one overlapping targets were then obtained as the related targets of GZFLW in the treatment of EMs (Fig. 3).

The data of the PPI network of those 101 targets were subsequently obtained in STRING Database. With confidence ≧0.7, there were 95 nodes and 751 edges in total. Taking two essential parameters of “degree” and “betweenness” as screening indexes, the topological analysis of targets mentioned above was performed. Targets more significant than or equal to the median are used as hub targets for GZFLW against EMs. The screened thresholds were degree ≧11.5 and betweenness ≧0.007, and the results were 25 hub nodes with 234 edges, including IL6, JUN, TNF, MAPK1, TP53, EGF, MAPK8, MMP9, VEGFA, AKT1, INS, FOS, ICAM1, PTGS2, CCL2, CCND1, EGFR, IL1B, MYC, IL10, PTEN, ESR1, PPARG, RELA and MMP2 (Table 2). When the 25 significant hub nodes and other 70 nodes were distributed with “degree” and “betweenness”, the network of 95 nodes was built as Fig. 4.

Based on the 25 key targets and related 29 active ingredients, we further established the network of Hub nodes—Compounds (Fig. 5). Among them, quercetin (MOL000098) is associated with 21 hub targets, and kaempferol (MOL000422) interrelates to 7 key targets. Besides, baicalein (MOL002714) is related to 6 key targets, while pachymic acid (MOL000289), ellagic acid (MOL001002), and taxifolin (MOL004576) act on 3 key targets respectively.

Molecular Docking

Thirteen hub targets with top degrees of GZFLW were identified with seven active compounds by AutoDock. There were 122 results (69.71%) of them had a moderate binding potential, which indicated that active ingredients of GZFLW could well bind to the targets for the treatment of EMs (Table 3). The ligands are mainly linked with corresponding proteins and critical amino acids around them in the form of hydrogen bonds (Fig. 6).

GO and KEGG pathway enrichment analysis

With the database of DAVID and Metascape, the enrichment analysis on 101 targets was performed and resulted in 116 GO items and 100 KEGG pathways.

GO enrichment analysis

After 116 items were sorted in descending order based on P-value, the first eight items of three parts, BP, MF, and CC, were selected (Fig. 7). In the aspect of BP, we mainly had: inflammatory response (GO:0006954), positive regulation of transcription from RNA polymerase II promoter (GO:0045944), transcription DNA-templated (GO:0006351), lipopolysaccharide-mediated signaling pathway (GO:0031663), positive regulation of cell division (GO:0051781), negative regulation of growth of symbiont in host (GO:0044130), response to toxic substance (GO:0009636) and positive regulation of cytokine secretion (GO:0050715); in the part of MF, we obtained heme binding (GO:0020037), sequence-specific DNA binding (GO:0043565), transcription factor activity sequence-specific DNA binding (GO:0003700), growth factor activity (GO:0008083), identical protein binding (GO:0042802), protein homodimerization activity (GO:0042803), cytokine activity (GO:0005125) and steroid binding (GO:0005496); in the aspect of CC, there were extracellular space (GO:0005615), nucleus (GO:00056340), cytosol (GO:0005829), membrane raft (GO:0045121), apical plasma (GO:0016324), external side of plasma membrane (GO:0009897), endoplasmic reticulum membrane (GO:0005789) and extrinsic component of external side of plasma membrane (GO:0031232). Based on the above three aspects, it is possible that the mechanism of GZFLW in treating EMs was the result of multi-pathway synergy.

KEGG pathway enrichment analysis

To further illustrate the potential mechanism of GZFLW in the treatment of EMs, we performed the KEGG pathway enrichment analysis on 101 targets. We selected the top 20 pathways based on the P-value, such as the AGE-RAGE signaling pathway(hsa04933), HIF-1 signaling pathway (hsa04066), PI3K-Akt signaling pathway (hsa04151), MAPK signaling pathway (hsa04010), EGFR tyrosine kinase inhibitor resistance (hsa01521) (Fig. 8). Then, we constructed the Target-Pathway Network to intuitively reveal the relationship between the hub targets and pathways (Fig. 9).

Primary endometrial stromal cells culture

The primary endometrial stromal cells were spindle-shaped or star-shaped with large and round nuclei, and proliferative fibrous tissue could be seen. After that, the immunofluorescence was used to identified the vimentin staining of cells in different groups (Fig. 10).

MTT colorimetric assay

In comparison with the control group, in addition to 0.5 mg/ml, the other different doses of treatment groups (1 mg/ml, 3 mg/ml, 5 mg/ml, 10 mg/ml, 20 mg/ml) had significant inhibitory effect on the proliferation of endometrial stromal cells (P < 0.05, 0.01 or 0.001; Fig. 11).

Western blotting

Western blot analysis revealed that compared to the control groups, drug treatment increased the protein levels of P53, BAX, and caspase3 (P < 0.05 or 0.01) in endometriotic lesions (Fig. 12).

Apoptosis analysis by flow cytometry

Compared with the control group, the percentage of apoptotic cells in the GZFLW group with the dose of 4 mg increased significantly (Fig. 13).

The major active ingredients of GZFLW

Based on the Hub nodes-Compound Network, we found a few hub compounds: quercetin and kaempferol contained in CM, taxifolin contained in CT, pachymic acid contained in PC, ellagic acid, baicalein, and paeoniflorin contained in RPR, campesterol, hederagenin contained in PK, β-sitosterol, contained in CM, RPR and PK. Among them, quercetin, kaempferol, taxifolin, and baicalein are flavonoids that have shown potent anti-EMs activity by effectively relieving symptoms and inhibiting levels of CA-125 (Signorile, Viceconte & Baldi, 2018). Quercetin and kaempferol can significantly reduce the ailing area of the endometrium via anti-proliferation and anti-inflammatory effects (Park et al., 2019) in EMs mice (Ilhan et al., 2020). Previous studies have demonstrated that phytosterols, including β-sitosterol and campesterol, contributed to the regression of EMs (Ilhan et al., 2019). Baicalin has also been reported to treat EMs by reducing the activity of endometrial stromal cell(ESC) (Jin, Huang & Zhu, 2017), which provided evidence for our results. In conclusion, we speculate that these above ingredients are the potential material basis of GZFLW in treating EMs.

It is mainly indicated that the pathways may be associated with treating EMs of GZFLW were regarding the immune response, apoptosis and proliferation, oxidative stress, and angiogenesis.

Inflammatory response

The differentially expressed genes are related to inflammation in EMs (Ahn et al., 2016), which suggest its role in the course of EMs. Prostaglandin-endoperoxide synthase (PTGS2), also known as COX-2, is related to the pain and infertility of EMs and the PTGS2/ PGE2 axis is considered as the critical target during EMs (Li et al., 2020). TNF, IL1- β are the other hub targets in this study. The immunoreactivity of which induces inflammation in the peritoneal fluid of EMs through activating NF-kappa B and mitogen-activated protein kinase (MAPKs) signaling pathways (Kralickova et al., 2018). To some extents, components such as hederagenin, campesterol, β-sitosterol, pachymic acid, quercetin, and taxifolin in GZFLW inhibit the expression of inflammatory factors to influence the whole inflammatory response. Previous studies have demonstrated that hederagenin inhibits the expression of inflammatory factors, including TNF, IL1, IL6, COX-2 via MAPKs and NF-kappa B pathway (Akhtar et al., 2019; Kim et al., 2017; Lu et al., 2015). Campesterol, β-sitosterol, known as sterols, can reduce the levels of TNF-α and IL of peritoneal fluid in the EMs mice (Ilhan et al., 2019). Flavonoids have a similar effect. Studies have experimentally proved that quercetin has an anti-inflammatory effect on the EMs (Park et al., 2019). Coincidentally, the docking result confirmed all the above, which also provides certain credibility for our results that GZFLW possibly worked against inflammatory in a multi-ingredient way.

Apoptosis and proliferation

The cell apoptosis’s inhibition and concomitant cell excessive proliferation play some important roles in the development of EMs. We think that the expression and regulation of the Jun and Fos genes are among the key factors in cell proliferation and apoptosis. It is found that high expression of c-Jun related to the proliferation of EMs cells via the JNK / c-Jun signaling pathway (Yu et al., 2018). TP53 is another critical tumor suppressor gene in our study; its regulative effect on the cell cycle plays a crucial role in obtaining proliferation activity for ESC of EMs (Hirakawa et al., 2016). It is investigated that the polymorphisms of TP53 may be involved in high risk of the generation of EMs (Hussain et al., 2018). Epidermal growth factor (EGF) also plays a vital role in regulating cell growth, proliferation, and differentiation. Previous studies have found that the EGF expression in patients with severe EMs significantly increased (Chatterjee et al., 2018). Furthermore, there are researches shown that COX-2 expression regulated by the PTGS2 gene can activate MAPKs, which then subsequently activate some transcription factors and protein kinases, ultimately promoting the proliferation of cells.

GZFLW’s regulative effect on cell proliferation and apoptosis in EMs is possibly achieved by regulating flavonoids on the expression of PTGS2 (Da Silva et al., 2020), which further regulates the COX-2/PGE2 axis (Takaoka et al., 2018). Besides, components such as flavonoids, quercetin in GZFLW affect the proliferation and apoptosis in EMs by inhibiting the ERK1/2, P38 MAPK, and AKT signaling pathway (Park et al., 2019). Also, baicalin can affect the activity of ESC by apoptosis inhibition via the NF-kappa B signaling pathway (Jin, Huang & Zhu, 2017). Furthermore, the direct or indirect inhibition of natural triterpenoids on proliferation is possibly another essential component of the regulation of cell apoptosis (Chen et al., 2019).

It is mentioning that all the primary active ingredients except paeoniflorin had a high docking score with EGF in our docking results. Besides, three active ingredients (hederagenin, phytosterol, and pachymic acid) had tight combinations with Jun, FOS, and PTGS2, and the other three active ingredients (pachymic acid, taxifolin, and hederagenin) had strong putative interaction with P53, which demonstrates that GZFLW is likely to affect on proliferation and apoptosis through a multi-ingredient synergistic way.

To verify the effect of GZFLW on EMs, we conducted in vitro experiments. The results of MTT and flow cytometry confirmed that GZFLW can effectively promote the apoptosis of ectopic stromal cells in EMs. In the treatment group, the expression levels of P53, caspase3 and BAX were significantly up-regulated, suggesting that the therapeutic effect of GZFLW may be related to the regulation of P53 pathway and apoptosis pathway.

Oxidative stress

More and more evidence has shown that the pathophysiological mechanism of EMs is complicated, and hypoxia plays an essential role in EMs damage (Wu, Hsiao & Tsai, 2019). The AGE/RAGE and HIF-1alpha signaling pathway were speculated as probable mechanisms of GZFLW against EMs from KEGG results. Previous studies have demonstrated that high expression of RAGE under hypoxia environment increased the accumulation of ROS in peritoneal fluid (Polak et al., 2018). Eventually, contribute to the cell invasion ability in EMs (Seguella et al., 2019). It is also found that high expression of HIF-1α in the serum of EMs aggravated the severity of dysmenorrheal (Zhang et al., 2018).

Fortunately, triterpenoids and flavonoids including hederagenin, β-sitosterol (Adebiyi et al., 2019; Ponnulakshmi et al., 2019), campesterol (Alvarez-Sala et al., 2018; Eom et al., 2017; Ho, Liu & Loke, 2016; Ogunlaja et al., 2016; Shamloo et al., 2017; Toiu et al., 2019; Weingartner et al., 2017), taxifolin, and pachymic acid were reported to be the natural antioxidants, which could reduce oxidative stress by inhibiting RAGE expression to a certain extent (Song et al., 2017) via MAPK, and PI3K/Akt pathways (Wen et al., 2018). The molecular docking results showed that triterpenoids and flavonoids have tight combinations with the hub targets AKT1, MAPK1, and MAPK8. However, further experiments are needed to confirm whether GZFLW worked against EMs by affecting AGE/RAGE and HIF-1α via MAPK and PI3K/Akt pathways.

Angiogenesis

Hypoxia and inflammatory, as an essential characteristic of EMs microenvironment, can stimulate the transcription of hypoxia-inducible factor HIF (Lin et al., 2018) and inflammatory-inducible factor COX-2 gene, consequently activate the angiogenesis through the PI3K/mTOR and MAPK pathways. HIF-1α/ VEGF expression in serum and endometrium has been reported to have a relationship with the stage of EMs and the severity of dysmenorrheal (Zhang et al., 2018). Studies confirmed that some plant inhibitors of COX-2 were able to treat EMs and inhibited ESC’s angiogenesis by regulating HIF-1α/VEGF in mice models (Li et al., 2020).

It is worth mentioning that triterpenoids (β-sitosterol, campesterol, hederagenin) and flavonoids (taxifolin) in GZFLW were detected in our study, they could have a modest effect on inhibiting VEGF in EMs, which have been confirmed by experiments and in line with our molecular docking results (Ilhan et al., 2019). Previous studies have also demonstrated that paeoniflorin could inhibit angiogenesis via the HIF-1 α/ VEGF pathway (Song et al., 2017; Zhou, Ding & Hardiman, 2018).