In our study, we designed a systematic pharmacological method to analyze the molecular mechanism of XFZYD in treating GBM (Figure 1). The key point of the method is to construct the weighted pathogenetic gene network and C-T network, then integrate the two networks through PPI data to construct a comprehensive C-T-P network, and then optimize the C-T-P network by using the novel node importance method to capture the key functional network. The dynamic programming algorithm was employed to predict CCG, finally, the CCG and their targets were used to infer the molecular mechanism of XFZYD in treating GBM. This strategy provided a methodological reference for formula optimization in TCM.
Construct Weighted Pathogenetic Gene network
The weighted pathogenetic gene network reflects the pathogenesis of glioma assisted by multiple genes. thus, the construction of pathogenic gene regulatory network is the key to provide evidence for glioma therapy. In this study, we use the pathogenetic genes (Figure S1) and comprehensive PPI network (Figure S2) to construct the weighted pathogenetic network. In this network, the highest weight gene is RRM2, following by TOP2A, MELK, NUSAP1, NDC80, DLGAP5, CHI3L1, IGF2BP3, METTL7B, and COL3A1. Most of these pathogenic genes were reported associated with cell apoptosis related pathways of GBM (Figure S3). Such as, Rikke et al. find that regulated RRM2 expression promoted tumorigenicity of GBM . Arivazhagan et al found that the expression of TOP2A was a prognostic biomarker for GBM patients treated with TMZ . Marie et al. confirmed that there exists a significant correlation between the expression of MELK and malignant grade of astrocytoma . Qian et al. demonstrated that NUSAP1 could be used as a new prognostic indicator and a underlying target for GBM . Clinical trials could be conducted to verify the efficacy of gene therapy by the down-regulating BUB1B (or NDC80) in GBM patients . These results show that the weighted pathogenetic genes network and weighted pathogenetic genes can reflect the pathogenesis of GBM, and also provide a reliable reference for the next step of constructing an intervention network (Figure S4).
Chemical components analysis
Chemical analysis is one of the main means to study the material basis and action mechanism of herbal medicine. Through literature search, we collected the chemical constituents and concentrations of various components in XFZYD. The detailed information was shown in Table 1. The results show that the chemical components of some herbs have high concentrations, suggesting that these components may have therapeutic effects and provide an extended experiment-aided chemical space for screening of active components and providing a more comprehensive reference for further analysis [3, 32, 33].
Among these document-confirmed components of XFZYD, paeoniforin [34-36], Gallic acid [37-39], Protocatechuic acid  and Chlorogenic acid [41-43] with more evidences that have an effect on GBM.
Active components of XFZYD
To fully understand the potential anti-GBM mechanism of XFZYD, we collected components of the 11 herbs in XFZYD from published databases, a total of 1365 components were collected. To assess the pharmacokinetic properties of the XFZYD formula, we used published ADME system to predict the pharmacological properties of BBB, OB, and DL. Only components with BBB＞-0.3, OB＞30%, and DL＞0.14 were predicted as active components. After screening, 117 active components were kept for further analysis. The detailed information of the active components can be found in Table 2.
C-T network construction
In order to observe the relation pattern of active components and their corresponding targets in XFZYD, Cytoscape was empolyed to construct C-T networks (Figure 2). The C-T network consists of 117 active components targeting to 855 genes through 4488 edges. Then the plug-in of Cytoscape, NetworkAnalyzer was used to analyze the topology parameters of th C-T network and showed that the average degree of targets and components in XFZYD were 5.25 and 39.40, respectively. The average number of components for per target is 5.25. This result indicates that XFZYD has multi-components feature in treating GBM. Among these components, gadelaidic acid (MOL004996, degree=123) has the highest number of targets, followed by icos-5-enoic acid (MOL004985, degree=115), cacetin (MOL001689, degree=107), 7-Methoxy-2-methyl isoflavone (MOL003896, degree=105), wogonin (MOL000173, degree=96), Glypallichalcone (MOL004835, degree=96), 1,3-dihydroxy-9-methoxy-6-benzofurano[3,2-c]chromenone (MOL004913, degree=90), Jaranol (MOL000239, degree=90), Mandenol (MOL001494, degree=86), and linolenic Acid (MOL000432, degree=84). Most of these components were reported associated with cell apoptosis related pathways of GBM. Such as cacetin (MOL001689, degree=107) showed lower IC50 values and presenting promising antit-glioma effects (Teles et al., 2015), wogonin (MOL000173, degree=96) induced human glioma cell apoptosis mediated reactive oxygen species (ROS) generation (Tsai et al., 2012), linolenic Acid (MOL000432, degree=84) was a selective activity against GBM cells compare normal brain(Das, 2007).
In the C-T network, the average number of targets for each component is 38.40. This result indicates that XFZYD has a multi-targets feature in treating GBM. In the top of 30 targets with larger weight include ESR2, ESR1, ABCG2, CYP19A1, PTN1, CA12, etc. Interestingly, the majority of these targets are related to pathogenesis or treatment of GBM. For example, ESR1 (degree=93) has GBM suppressor function and ERβ isoform switching induces GBM progression (Liu et al., 2018), ABCG2 (degree=81)
was an important indicator of patient prognosis, and ABCG2 inhibitor may be a potential treatment strategy for GBM (Emery et al., 2017)(Liu and Sareddy et al., 2018), the CA12 (degree=26) inhibitor prolongs the survival of TMZ-treated GBM-bearing mice . Thus, these components could be considered as curative elements in treating GBM.
Key Functional Network Identification
The multi-component and multi-target characteristics of TCM determine that it exists network intervention in the process of treating complex diseases, and the active components in the formula can transmit the intervention effect to pathogenic genes through PPI networks by influencing their targets, thus playing a therapeutic role on diseases. Components-targets-pathogenic genes form a complex intervention network, in which some components play a crucial role and some components play an auxiliary role. How to find the core component groups and key functional networks from this network is the foundation to understand the treatment of complex diseases of formula in TCM and is also the foundation for the secondary development of TCM formula. PageRank (PRK) algorithm usually used in the Google search engine to measure the importance of web pages, is also widely used in biomedical researches. Syed et al. used PRK centrality to consistently find convincing eloquent brain areas (Perry et al., 2019). Hao Zhang et al. predicted human and plant target genes using RNAhybrid by PRK algorithm (Zhang et al., 2016). In this paper, we combined C-T network, weighted pathogenetic gene network and PPI network to form a comprehensive network, in this network, we use PRK value to describe the importance and influence of nodes, the nodes with higher PRK value than the average and median PRK value were defined as A functional network and M functional network, respectively. In the further analysis of A functional network and M functional network, we found that A functional network includes 37 components and 728 targets, M functional network includes 33 components and 688 targets. Then we performed pathway analysis by using targets of A functional network, M functional network and C-T network respectively. Among them, the number of targets in A functional network enriched pathways was 131 (p < 0.05), the number of targets in M functional network enriched pathways was 125 (p < 0.05), the number of genes in C-T network enriched pathways was 154 (p < 0.05). Compared with enriched pathways of genes in the C-T network, the coverage rates of the two strategies were 77.92% and 74.68%, respectively (Figure 3). According to this result, we used the A functional network as key functional network for further analysis.
CCG Prediction and validation
The key functional network contains a large number of chemical components, action targets and interactions between targets, which form a complex response network. There are potentially effective component groups and modes of action in this network. How to detect these effective component groups and action patterns is key to understand the treatment of glioma with XFZYD. Here we design CCR model to obtain CCG by heuristic optimization of key functional networks using reasonable global metrics. According to the results of CCR model, the top 5 components including gadelaidic acid (MOL004996), icos-5-enoic acid (MOL004985), acacetin (MOL001689), 7-Methoxy-2-methyl isoflavone (MOL003896), and Glypallichalcone (MOL004835) contribute to 55.63% target coverage of genes in key functional network. For further analysis, 21 components can contribute to 90.66% targets coverage of genes in key functional network were selected as CCG (Figure 4 and Table 3). Higher targets coverage of genes in key functional network proved that the CCG may play a leading role and produced combination actions in treating GBM.
To analysis of XFZYD in the treatment of GBM at the functional level, we used CCG targets for pathway analysis, the number of CCG targets enriched pathways was 121 (p < 0.05), and the number of pathogenic genes enriched pathways was 26 (p < 0.05). We defined the shared enrichment pathways of pathogenic genes and targets in key functional network as the reference pathways. The CCG targets enriched pathways account for 78.95% of the reference pathways (Figure 5A). This results confirmed the accuracy and reliable of our CCG selection model. For futher analysis, we found that these main targets of CCG were frequently involved in PI3K/Akt signaling pathway (hsa04151), Toll-like receptor signaling pathway (hsa04620), Proteoglycans in cancer (hsa05205), Cell cycle (hsa04110), Cellular senescence (hsa04218), and Platelet activation (hsa04611), etc (Figure 5B). For example, PI3K/AKT pathway plays a key role in the occurrence and development of GBM . AKT, as a regulatory molecule and a potential drug target, has been widely studied for its carcinogenicity [46, 47]. AKT activation promotes tumor progression by affecting tumor proliferation and growth . Up to now, three isomers of AKT (AKT1, AKT2 and AKT3) have been discovered. According to reports, AKT3 can significantly activate DNA repair to resist radiotherapy and chemotherapy in GBM patients (Turner et al. (2015). Results show that our combined optimized strategy of key functional network detection and CCR model is accurate and reliable, and our predicted core component group plays an critical role in treating GBM.
GO enrichment analysis of CCG targets
To further analyze the combined effects of XFZYD, all targets of CCG were enriched by GO enrichment analysis (Figure 6B). We defined the shared enriched GO terms of key functional network and pathogenetic gene as the therapeutic response GO terms. The enriched GO terms of CCG targets account for 73.68% of the therapeutic response GO terms (Figure 6A).
GO enrichment analysis indicated that the targets regulated by CCG of XFZYD were mainly enriched in oxidative stress and lipid metabolism. For example, the GO terms of oxidative are response to oxidative stress (GO:0006979, CYP1B1, LCN2, PTGS1, FABP1, PTGS2, CDK1, etc.), decreased oxygen levels (GO:0036293, CYP1A1, PRKCE, RORA, ADORA1, PPARA, PPARD, TERT, etc.), hypoxia (GO:0001666, CYP1A1, PRKCE, RORA, ADORA1, PPARA, PPARD, TERT, etc.), and cellular response to oxidative stress (GO:0034599, CYP1B1, LCN2, FABP1, CDK1, MCL1, NOX4, MMP2, etc.). Important genes involved in oxidative stress in these GO terms include PPARG, PPARA, PPARD, TERT, CDK1, etc. These genes are involved in the oxidative stress of GBM by fostering proliferation and affecting survival and migration [49-51]. The GO terms of lipid metabolism are regulation of lipid metabolic process (GO:0019216, CYP1A1, NCOA1, PRKCE, FABP3, NR1H3, RORC, HMGCR, CYP51A1, etc.), lipid transport (GO:0006869, ABCB1, CYP19A1, ABCC1, NCOA1, SLCO2B1, PLA2G2E, FABP4, FABP3, etc.), lipid localization (GO:0010876, ABCB1, CYP19A1, ABCC1, NCOA1, SLCO2B1, PLA2G2, etc.), response to lipopolysaccharide (GO:0032496, CYP1A1, PRKCA, LCN2, MAOB, PRKCE, NR1H3, PPARD, etc.) and lipid catabolic process (GO:0016042, FABP7, CYP1B1, CYP19A1, FAAH, PLA2G2E, SRD5A1, PRKCE, etc.). Important genes involved in lipid metabolism in these GO terms include ABCB1, ABCC1, FABP4. These genes are involved in the lipid metabolism of GBM by influencing cell growth, proliferation, angiogenesis and invasion [52, 53]. GO analysis confirmed that XFZYD treats GBM by regulating of oxidative stress and lipid metabolism.
Pathway analysis to explore the potential therapeutic mechanisms of XFZYD
In order to investigate the potential therapeutic mechanisms of XFZYD on GBM, We analyzed the intersection of enrichment pathways of CCG target gene and GBM pathogenic genes, totally get 15 shared pathways. Among the 15 common pathways, remove the pathways associated with other diseases, the PI3K-Akt signaling pathway (hsa04151) and Toll-like receptor signaling pathway (hsa04620) were the literature proved GBM-related pathways. For example, PI3K-Akt signaling pathway plays an important role in cell cycle, apoptosis, and cell proliferation in GBM [54, 55].
In order to analysis the synergetic mechanism of XFZYD in treating of GBM, we constructed a comprehensive signaling pathway by combining two principal molecular pathways. As shown in Figure 7, we defined the first three columns as the upstream position of the pathway and the remaining columns as the pathway downstream position. Among them, targeting to the PI3K-Akt signaling pathway is the main approach for XFZYT to treat GBM. XFZYD regulates 7 targets located the upstream, such as Syk, RTK, TLR2 / 4, and 32 targets located PI3K-Akt signaling pathway downstream, such as MEK, AKT, and mTOR. Downstream targets account for 82.05%. XFZYD may activate downstream of AKT proteins through upstream RTK, leading to cascade amplification of downstream GSK3 and CREB, which was closely related to GBM cell proliferation and protein synthesis .
Toll-like receptor signaling pathway is also a critical way for XFZYD to treat GBM. The target of XFZYD regulation is located the pathway downstream. Such as, 11 targets of XFZYD, including AKT are located downstream of the Toll-like receptor signaling pathway. XFZYD can affect upstream TLR2/4 and Rac1 and then activate downstream AKT, thereby further affecting a series of malignant biology with GBM behavior-related proteins, such as STAT1, NFκB, and IL-6, thus play an important role in the treatment of GBM.
Experimental Validation in Vitro
The effects of baicalein and wogonin in CCG with different concentrations on the viabilities of human LN229 glioma cells were detected by CCK8 assay. After 72h of incubation, the cell viabilities of LN229 cells were 79.63%, 67.68%, 50.38%, 44.61%, 36.28% and 88.61%, 78.91%, 66.88%, 53.36%, 44.53% after exposure to baicalein and wogonin at concentrations of 40, 60, 80, and 100 μM, respectively (Figure 8A). The above results demonstrated that baicalein and wogonin remarkably suppressed cell proliferation in LN229 cells.
According to the results of cell proliferation effect of baicalein and wogonin, we found 60 μM baicalein and 80 μM wogonin could markedly inhibit the proliferation of LN229 cells. Therefore, 60 μM baicalein and 80 μM wogonin were used to evaluate the effect on apoptosis in LN229 cells for 24 h. Compared with the control group, Necrosis and late apoptotic cells of upper right quadrant was increasing (65.1% vs 0.59%, 52.2% vs 0.59%, respectively) with baicalein and wogonin treatment (Figure 8B). Further, we found that treatment with baicalein and wogonin could significantly reduced Bcl-2 protein expression (Figure 8C). Our data suggested that the baicalein and wogonin induced apoptosis phenomenon play a important role in inhibiting the development of GBM.