A System Pharmacology Model for Decoding the Synergistic Mechanisms of Compound Kushen Injection in Treating Breast Cancer.

Breast cancer (BC) is one of the most common malignant tumors among women worldwide and can be treated using various methods; however, side effects of these treatments cannot be ignored. Increasing evidence indicates that compound kushen injection (CKI) can be used to treat BC. However, traditional Chinese medicine (TCM) is characterized by "multi-components" and "multi-targets", which make it challenging to clarify the potential therapeutic mechanisms of CKI on BC. Herein, we designed a novel system pharmacology strategy using differentially expressed gene analysis, pharmacokinetics synthesis screening, target identification, network analysis, and docking validation to construct the synergy contribution degree (SCD) and therapeutic response index (TRI) model to capture the critical components responding to synergistic mechanisms of CKI in BC. Through our designed mathematical models, we defined 24 components as a high contribution group of synergistic components (HCGSC) from 113 potentially active components of CKI based on ADME parameters. Pathway enrichment analysis of HCGSC targets indicated that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis could synergistically target the PI3K-Akt signaling pathway and the cAMP signaling pathway to treat BC. Additionally, TRI analysis showed that the average affinity of HCGSC and targets involved in the key pathways reached -6.47 kcal/mmol, while in vitro experiments proved that two of the three high TRI-scored components in the HCGSC showed significant inhibitory effects on breast cancer cell proliferation and migration. These results demonstrate the accuracy and reliability of the proposed strategy.

Introduction

Breast cancer (BC) is the most common cancer in women worldwide and is responsible for the second highest death rate among female patients with cancer (Siegel et al., 2020). Many options are available to treat breast cancer, such as surgical treatment (Zheng et al., 2020), radiation therapy (Balaji et al., 2016), neoadjuvant endocrine therapy (Yao et al., 2019), neoadjuvant chemotherapy (Vaidya et al., 2018), and anti-HER2 therapy (Yao and Fu, 2018). However, side effects of these treatments are often observed. For example, radiation therapy can cause radiation-induced fibrosis (Straub et al., 2015) and radiodermatitis (Singh et al., 2016). Fatigue, pain, and systemic side effects can occur as a result of surgical treatment (Enien et al., 2018). Adjuvant endocrine therapy can cause vasomotor symptoms and musculoskeletal and vulvovaginal symptoms (Condorelli and Vaz-Luis, 2018). These side effects mainly influence the organs and blood system, which aggravates the psychological burden on the patients. Side effect symptoms include loss of appetite, vomiting, diarrhea, nausea, and ulcers (Meng et al., 2017; Recht, 2017; Cochran et al., 2019). In recent years, traditional Chinese medicine (TCM) has become increasingly popular in the treatment of BC. Ganoderma lucidum from Ganoderma suppresses the proliferation and migration of breast cancer by inhibiting Wnt/β-catenin signaling (Zhang, 2017), while puerarin from Radix puerariae inhibits cell migration, invasion, and adhesion of LPS-induced BC by blocking NF-κB and Erk pathways (Liu et al., 2017). In addition to these herbal components, many formulas and proprietary Chinese medicines are widely used in the treatment of BC, such as compound kushen injection (CKI) (Liu et al., 2020), Shugan Jianpi decoction (Jingyuan et al., 2021), and Fangjihuangqi decoction (Guo et al., 2020). Based on these formulas and proprietary Chinese medicines, CKI is a commonly used antitumor treatment in clinical practice.

CKI is a TCM formulae extract from kushen (Radix Sophorae flavescentis) and baituling (Rhizoma Heterosmilacis) at a ratio of 7:3, in which exist hundreds of components including alkaloids and flavonoids (Cui et al., 2020). As approved by the China Food and Drug Administration, CKI can be employed for cancer treatment (Guo et al., 2015) and has several pharmacological functions, including anticancer properties, hemostasis, and immunity enhancement (Wang et al., 2015). CKI is widely being used in the treatment of cancers, such as breast cancer (Xu et al., 2011a; Qu et al., 2016; Cui et al., 2019; Nourmohammadi et al., 2019; Cui et al., 2020), acute myeloid leukemia (Jin et al., 2018), and hepatocellular carcinoma (Gao et al., 2018; Yang et al., 2021). Additionally, the function of CKI in treating breast cancer was also proved in in vitro and in vivo experiments. CKI could reduce the tumor formation rates and tumor volume by the downregulated Wnt/b-catenin pathway in the MCF-7 SP xenograft model (Xu et al., 2011a). Besides, CKI could impair the migration and invasiveness of MDA-MB-231 cell lines (Nourmohammadi et al., 2019) and influence the cell cycle of MDA-MB-231 cell lines in promoting cell death by increasing the proportion of cells in the G1 phase and decreasing the cells in the S and G2/M phase (Cui et al., 2019; Cui et al., 2020), while it could inhibit proliferation and induce cell apoptosis of MCF-7 cell lines (Qu et al., 2016). These in vitro and in vivo experiments had proved that CKI has anticancer pharmacological function, especially in breast cancer. In CKI, matrine, oxymatrine, and sophocarpine are the main components of Radix Sophorae flavescentis, which display antitumor, anti-inflammatory, and antiviral properties, as well as cardiovascular protective abilities (Wang et al., 2015). Rhizoma Heterosmilacis may play a role in promoting oxidative stress–induced apoptosis (Liu et al., 2017), reducing oxidative stress (Hong et al., 2014), deoxidation, and dampness relief (Liang et al., 2019). Increasing evidence has shown that CKI can be used to halt cancer migration (Nourmohammadi et al., 2019), reduce anticancer drug resistance (Hy et al., 2014), increase cancer cell apoptosis (Qu et al., 2016), suppress the cancer cell cycle (Cui et al., 2019), and inhibit cancer progression (Wang et al., 2019). Matrine upregulates Bax and downregulates Bcl-2 to inhibit proliferation and increase apoptosis of breast cancer cells (Li et al., 2015) and further downregulates the canonical pathway to suppress human breast cancer stem-like cells (Xu et al., 2011b). Oxymatrine exerts its effects by blocking the cell cycle, initiating apoptosis (Binggang et al., 2002) and suppressing the epithelial–mesenchymal transition (Jiaqin et al., 2018). These sporadic reports suggest that different components of CKI play a role in the treatment of BC, but there is still a lack of systematic and overall research on the synergistic mechanism of different components of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis in CKI.

Since TCM is characterized by “multi-components” and “multi-targets”, it is difficult to reveal the associated mechanisms between herbs, components, genes, and disease through traditional experimental methods. As an effective tool to elucidate the synergistic and potential mechanisms of the networks between component–target and target–disease, systemic pharmacology provides a new perspective on the therapeutic mechanisms of TCM in the treatment of BC at the systematic level.

To explore the therapeutic mechanism of CKI in treating BC, a novel system pharmacology strategy (Figure 1), which integrates differentially expressed gene analysis (DEGs), pharmacokinetics synthesis screening, target identification, synergy contribution degree (SCD), network analysis, and therapeutic response index (TRI) calculation, was developed. The Cancer Genome Atlas (TCGA) database was used to explore the DEGs in BC. Subsequently, previously proposed absorption, distribution, metabolism, and excretion (ADME) screening models were employed to select potential active components; online tools were used to predict the targets of these potential active components, and a component–target (C-T) network was used for further network analysis. Network analysis of potential components combined with SCD was used to simulate the treatment effects of each component of CKI in BC and to explore a high contribution group of synergistic components (HCGSC). The virtual docking-aided TRI model was designed to determine whether highly reliable components may have a potential synergistic mechanism. Thereafter, GO enrichment and KEGG pathway enrichment analysis of BC, CKI, and highly reliable components in HCGSC were discussed to decode the potential synergistic mechanism analysis of CKI in the treatment of BC based on HCGSC. Ultimately, experimental validation was used to confirm the effect of high TRI-scored components in the HCGSC and evaluate the reliability of our model. We hope that these results will provide a strategy to reveal the therapeutic mechanism of TCM at the molecular level.

FIGURE 1

Workflow of the system pharmacology approach.

Materials and Methods

Chemical Components Collected From Database

All CKI components were collected from the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (http://www.tcmspw.com/tcmsp.php). TCMSP includes information from DrugBank, HIT, TTD, and PharmGKB. The pharmacokinetic properties of TCMSP include molecular weight (MW), AlogP, number of acceptor atoms for H-bonds (nHacc), number of donor atoms for H-bonds (nHdon), Caco-2 permeability (Caco-2), oral bioavailability (OB), drug-likeness (DL), blood–brain barrier (BBB), FASA-, and half-life (HL).

ADME Screening

In the process of modern drug development, many drug candidates fail during development because of inadequate ADME properties (Rojas-Aguirre and Medina-Franco, 2014). Therefore, evaluating the ADME of drug components in the early stages has become an essential process. Components with better pharmacokinetic properties can be obtained by ADME screening, and the potential drug–drug interactions can be minimized (Wang et al., 2018). Two ADME-related models were employed in the present study, namely Caco-2 permeability and drug-likeness, to screen the potential active components of CKI (Supplementary Figure S1).

As a model system for intestinal epithelial permeability, the human colon carcinoma cell line (Caco-2) (Ij et al., 1989) is currently considered the gold standard (Jacobsen et al., 2020). The transport rates of components (nm/s) in Caco-2 monolayers represent intestinal epithelial permeability in TCMSP (Ru et al., 2014). Components with Caco-2 < −0.4 are not permeable; thus, Caco-2 > −0.4 were selected as candidate components.

Drug-likeness (DL) may be defined as a complex balance of various molecular properties and structural features that determine whether a particular molecule is similar to known drugs. The drug similarity index of the new compound was calculated using the Tanimoto coefficient, which is defined as

In this equation, A represents the molecular descriptor of herbal components, and B is the average molecular property of all components in DrugBank (Tao et al., 2013). Based on the data from TCMSP, the average value of Radix Sophorae Flavescentis was 0.40, while that of Rhizoma Heterosmilacis was 0.25. Combined with other literature search criteria, we defined DL ≥ 0.18, as the screening criterion of DL.

Target Identification

Three commonly used online tools were employed to identify the targets of active components in CKI, namely, similarity ensemble approach (SEA) (Keiser et al., 2007), HitPick (Liu et al., 2013), and SwissTargetPrediction (Gfeller et al., 2014). Open Babel 3.0.0 (O Boyle et al., 2011) was used to convert the SDF format of all CKI potential active components into the SMILES format. Potential active components in the SMILES format were imported to the SEA, HitPick, and SwissTargetPrediction to predict targets of potential active components. The homologous genes of other species provided by online tools were included in the discussion.

Network Construction and Analysis

The C-T network was used as a frame to uncover the relationship between active components and targets. Cytoscape 3.8.0 (Shannon et al., 2003), an open-source software platform, was employed to visualize the networks.

GO Enrichment and KEGG Pathway Enrichment Analysis

To analyze the main function of targets, clusterProfiler (Yu et al., 2012) in the Bioconductor package (https://bioconductor.org/) based on the R language was used for GO-BP enrichment analysis and KEGG pathway analysis. FDR-adjusted p values were set at 0.05, as the cut-off criterion.

Synergy Contribution Degree Calculation

The SCD represents the contribution of a potential active component in the C-T network and its effectiveness in the treatment of BC. To evaluate the effect of CKI on the treatment of BC, we built a mathematical model to calculate the SCD of each active component in CKI:

In this equation, i is the number of components, and j is the number of targets. C represents the betweenness centrality of each component. T represents the betweenness centrality count of all targets of each component. Ce represents the edge count of each component. Te represents the edge count of all targets of each component. C represents the betweenness centrality of each component only in Rhizoma Heterosmilacis, and C represents the betweenness centrality of each component only in Radix Sophorae Flavescentis. If a component consists only of Radix Sophorae Flavescentis or Rhizoma Heterosmilacis, should be 1. If is zero, we assign 1e-6 to it. d represents the dose of each component in CKI, which was extracted using HPLC (Supplementary Table S1). NORM (di) is the index of min-max normalization to the dose of each component. The betweenness centrality of each component and target was calculated using Cytoscape version 3.8.0. based on the C-T network of the potential active components.

Breast Cancer Differentially Expressed Genes Analysis

The TCGA database (https://portal.gdc.cancer.gov/) was used to explore the DEGs in BC. We used TCGAbiolinks (Colaprico et al., 2016) in the Bioconductor packages to access the HTSeq-Counts number of normal and BC samples from the TCGA-BRCA project of the TCGA program. DESeq (Anders and Huber, 2010) and limma (Ritchie et al., 2015) in the Bioconductor package were used to analyze the DEGs in BC. We set |log2FoldChange|≥2 and FDR adjusted the p value < 0.05, as the cut-off criterion. The results are presented and plotted using the gplots (Warnes et al., 2020) package.

Molecular Docking Analysis of HCGSC in the Treatment of BC

The components of HCGSC were collected from the ZINC (Sterling and Irwin, 2015) and PubChem (https://pubchem.ncbi.nlm.nih.gov) databases in the MOL2 format. Human proteins were collected from the Protein Data Bank (PDB) (http://www.rcsb.org). AutoDock Tools (ADT) (Sanner, 1999) was used to pretreat the components and proteins. The pocket of the protein was automatically extracted from the space of the protein. AutoDock Vina (Trott and Olson, 2009) was used for docking. The seed of docking was 10,000, the energy range was 4, and the exhaustiveness was 4. The affinity (kcal/mol) index of each component–protein pair was used to estimate the docking results. The results were obtained using PyMOL (Schr Odinger, 2015).

Therapeutic Response Index Model Calculation

We developed a mathematical model to calculate the therapeutic response index (TRI) of each component in a highly reliable docking relationship (HRDR) between HCGSC and its targets.

In the above model, i is the number of components, j is the number of topology parameters, n is the collection of components, and m is the collection of components and topology parameters. CDI represents the docking index of HRDR. GI represents the number of genes targeted by each component. AI represents the average affinity of all component–protein binding relationships of each component. NI represents the C-T network topology parameters of HCGSC and their targets, including degree, betweenness centrality, closeness centrality, average shortest path length, eccentricity, radiality, neighborhood connectivity, and topological coefficient. SCD is the calculation result for each component from method 2.6. TRI represents the mechanism of CKI in the treatment of BC, which is the sum of normalized CDI and normalized SCD. NORM is a normalized equation.

Materials

High-glucose Dulbecco’s modified Eagle medium (DMEM; 4.5 g/L), fetal bovine serum (FBS), 0.25% trypsin, and 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) were purchased from Shanghai Sangon Biotechnology Co., Ltd. (Shanghai, China). The transwell invasion chamber was purchased from Corning costar company (Andover, MA, United States). Matrigel was purchased from BD (Biosciences, Bedford, MA, United States). Luteolin and trifolirhizin (≥98% purity by HPLC) were obtained from Nanjing Jingzu Biotech Co., Ltd. (Nanjing, China).

Cell Culture

The human breast cancer cell line MCF-7 was purchased from Procell Life Science & Technology Co., Ltd. Cells were cultured in an incubator at 37°C in DMEM containing 10% FBS. When the cells reached 80% confluency, they were exposed to different concentrations of luteolin and trifolirhizin (1, 2, and 4 mg/ml) for 24 h.

Cell Viability Assay

MCF-7 cells (1 × 104 cells/well) were seeded in 96-well plates and treated with 0, 5, 10, 20, 40, 80, and 100 μM luteolin and trifolirhizin for 24 h. MTT was added to a 96-well plate for 4 h, and the culture supernatant was removed. Finally, DMSO was used to dissolve the purple crystals. Absorbance at 570 nm was measured using a plate reader.

Transwell Assay

Referring to our previous method (Gao et al., 2018), luteolin (40 and 80 μM) and trifolirhizin (20 and 40 μM) were added to the lower compartments. The migrating cells were observed under a microscope.

Statistical Analysis

The R package ggpubr (Kassambara, 2020) was used to compare the molecular properties of all components in Rhizoma Heterosmilacis and Radix Sophorae Flavescentis. Data were analyzed using the Student’s t-test. Differences were considered statistically significant at p < 0.05.

Results

Based on the systematic pharmacological model, the mechanism of CKI in BC treatment was clarified and validated. First, DEGs that could be used as potential pathogenic genes were determined based on the data from the TCGA-BRCA project. Second, the CKI components were collected from the TCMSP database, and the ADME method was used to screen active components in CKI. Third, the active components and their targets, which were predicted using three online tools, were used to construct the C-T network. Fourth, a network analysis of potential components combined with SCD was used to explore HCGSC. Thereafter, GO enrichment and KEGG pathway enrichment analysis of BC, CKI, and highly reliable components in HCGSC were discussed to decode the potential synergistic mechanism of CKI in the treatment of BC. Subsequently, the virtual docking-aided TRI model was designed to determine whether highly reliable components may have a synergistic mechanism. Ultimately, experimental validation was used to confirm the effect of high TRI-scored components in HCGSC and evaluate the reliability of our model.

Differentially Expressed Genes Analysis of BC

To further analyze the mechanism of CKI in the treatment of BC, 113 normal breast samples and 1,102 breast cancer samples were extracted from the TCGA-BRCA project. A total of 979 DEGs were upregulated, and 990 were downregulated in cancer (Supplementary Table S2). The DEGs were used to determine the gene expression patterns of normal and tumor patients (Supplementary Figure S2). The results showed that the DEG pattern distinguished between diseased and normal states. From the expression pattern, we identified that the top 10 upregulated genes were FTHL17, CSAG1, MUC2, COX7B2, CGA, CSAG4, CST4, MAGEA12, MAGEA1, and ACTL8. Among these genes, MUC2 influences proliferation, apoptosis, and metastasis of breast cancer cells (Astashchanka et al., 2019), while the upregulation of MAGEA in patients revealed a higher risk of recurrence (Otte et al., 2001). The top 10 downregulated genes were LEP, GLYAT, AC087482.1, APOB, TRHDE-AS1, AQP7P2, FP325317.1, PLIN1, CA4, and AL845331.2. LEP which inhibit apoptosis of breast cancer cells by coding leptin, while PLIN1 inhibits invasion, migration, and proliferation of cells (Zhou et al., 2016; Crean-Tate and Reizes, 2018). The literature reports proved that these genes are related to the development of BC, which indicates that the DEGs were more likely to be pathogenic genes.

Component Comparisons in Rhizoma Heterosmilacis and Radix Sophorae Flavescentis

From the TCMSP database, 187 compounds were retrieved from Rhizoma Heterosmilacis (74) and Radix Sophorae Flavescentis (113). Detailed information on the components of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis is provided in Supplementary Table S3. To further describe the differences between Rhizoma Heterosmilacis and Radix Sophorae Flavescentis, we compared nine properties of the two components, including drug likeness (DL), Caco-2, molecular weight (MW), AlogP, TPSA, FASA-, RBN, nHDon, and nHAcc, and drew boxplots using R packages ggplot2 (Wickham, 2016), ggpubr (Kassambara, 2020), and ggsci (Xiao, 2018). The results (Figure 2) showed most of chemical composition and properties between Rhizoma Heterosmilacis and Radix Sophorae Flavescentis had no statistical difference except for DL (p = 0.00013) and RBN (p = 0.015). These results indicate that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis have similar chemical properties.

FIGURE 2

Boxplot of nine properties of the two components. The box in blue represents the value distribution of Radix Sophorae Flavescentis, and the box in red represents the value distribution of Rhizoma Heterosmilacis.

Potential Active Components in Rhizoma Heterosmilacis and Radix Sophorae Flavescentis

Traditional Chinese medicine uses a combination of herbs, each of which usually contains hundreds of components, and only a few of these components possess satisfactory pharmacodynamic and pharmacokinetic properties. In this study, two ADME-related properties, DL and Caco-2, were used to screen for active components. Through ADME screening, 111 active components (Radix Sophorae Flavescentis 90 and Rhizoma Heterosmilacis 21) were selected from the 187 components of CKI. Additionally, some components, such as N-methylcytisine and trifolirhizin, did not pass the ADME screening, but were frequently reported in previous studies, so we added them to the list of active components (Xiumei and Cen, 2004; Juan et al., 2007; Yue, 2012; Liang et al., 2013). Therefore, 113 active components (Radix Sophorae Flavescentis 92 and Rhizoma Heterosmilacis 21) were included for further analysis. Detailed information is provided in Supplementary Table S4.

Through ADME screening and literature review, 92 potential components containing ideal pharmacokinetic profiles were selected from Radix Sophorae Flavescentis. For example, as one of the major components of CKI, sophoridine (KS18, DL = 0.25, Caco-2 = 1.13) displays antitumor, anti-allergic, anti-inflammatory, anti-arrhythmia, and antiviral properties, and affects the central nervous system (Han et al., 2009; Cuiping, 2010; Huang et al., 2014; Hu et al., 2016). Anagyrine (KS34, DL = 0.24, Caco-2 = 1.16) and (+)-lupanine (KS83, DL = 0.24, Caco-2 = 1.16) inhibit the proliferation and induction of apoptosis in human cervical cancer cells (Merghoub et al., 2011). Oxysophocarpine (KS89, DL = 0.29, Caco-2 = 1.042) suppressed cell proliferation, migration, invasion, and angiogenesis, induced cell cycle arrest, and enhanced apoptosis in oral squamous cell carcinoma (Liu et al., 2018). IDO1 expression can be downregulated by kushenol E (KS76, DL = 0.59, Caco-2 = 0.58) and kushenol F (KS74, DL = 0.61, Caco-2 = 0.45) to inhibit tumor proliferation (Kwon et al., 2019) and induce apoptosis (Kwon et al., 2020).

Among the 74 components in Rhizoma Heterosmilacis, 21 components met the screening criteria. For instance, diosgenin (BTL5, DL = 0.81, Caco-2 = 0.82), a steroid compound, has been shown to promote apoptosis and anticancer effects (Romero-Hernández et al., 2015; Jiang et al., 2016). Sitogluside (BTL1, DL = 0.62, Caco-2 = -0.14), beta-sitosterol (BTL2, DL = 0.75, Caco-2 = 1.32), sitosterol (BTL3, DL = 0.75, Caco-2 = 1.32), and stigmasterol (BTL4, DL = 0.75, Caco-2 = 1.45) can inhibit proliferation (Awad et al., 2000; Awad et al., 2003) and promote apoptosis (Awad et al., 2003), cell cycle arrest, and sphingomyelin cycle activation (Jia et al., 2018). Additionally, quercetin (BTL1, DL = 0.28, Caco-2 = 0.05) could be used to inhibit the proliferation of tumor cell lines (Ren et al., 2015) and induce apoptosis of tumor cells (Jia et al., 2018). Taxifolin (BTL13, DL = 0.27, Caco-2 = -0.23) inhibited proliferation, migration, and invasion of breast cancer cells by promoting EMT through β-catenin signal transduction (Li et al., 2019).

Target Prediction and Analysis of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis

To evaluate the effect of the active components of CKI, we used SEA, HitPick, and SwissTargetPrediction tools to predict the targets of 113 active components of CKI. Finally, 780 targets were obtained from these tools (Supplementary Table S5). To explore the mechanism of CKI in the treatment of BC, 113 active components and 780 targets were used to construct the C-T network (Supplementary Figure S3). The C-T network results showed that 252 of all 780 targets were overlapped by the target list of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis, and 320 targets were unique to Radix Sophorae Flavescentis, while 208 targets were unique to Rhizoma Heterosmilacis. A total of 4,592 component–target associations between 113 active components and 780 targets were contained in the C-T network. The average number of targets per component is 40.64, and the mean number of components per target is 5.89, which shows that CKI has multi-component and multi-target characteristics in treating BC.

Among all components, six components displayed a number of degrees above 120, namely palmitone (BTL15, degree = 149), quercetin (BTL1, degree = 144), quercetin (KS3, degree = 144), norartocarpetin (KS77, degree = 130), apigenin (KS2, degree = 129), and luteolin (KS1, degree = 125), all of which exhibit anticancer functions (Imran et al., 2019), and suppress migration and EMT (Cao et al., 2020). These six high-degree components only accounted for 5.3% of all components but covered 58.3% of all common targets. They all targeted common targets with high degrees, including ESR1 (degree = 65), MAPT (degree = 65), and ESR2 (degree = 63). It is worth noting that the top 25 targets with the highest degree were all dropped in the common targets list between Rhizoma Heterosmilacis and Radix Sophorae Flavescentis, most of which were related to the pathogenesis or treatment of BC. As the coding genes of the estrogen receptor, methylation of the ESR1 (degree = 65) promoter may be associated with shorter survival time and increased risk of drug resistance to anti-endocrine therapy (Kirn et al., 2018) and is commonly targeted by oxysophocarpine (KS89), kushenol E (KS76), and kushenol F (KS74) of Radix Sophorae Flavescentis and Stigmasterol (BTL5), quercetin (BTL1), and taxifolin (BTL13) of Rhizoma Heterosmilacis. Additionally, MAPT (degree = 65) is correlated with microtubule assembly and stabilization (Weingarten et al., 1975), which can promote bicalutamide resistance and is associated with survival in prostate cancer (Sekino et al., 2020). Our data analysis showed that (+)-Lupanine (KS76) and leachianone,g (KS74) from Radix Sophorae Flavescentis and quercetin (BTL1) from Rhizoma Heterosmilacis could target MAPT. We noted that several components of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis can both activate common targets to play their role in the treatment of BC. Thus, these results suggest that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis may act synergistically to treat BC via the common targets based on the “multi-component” and “multi-target” features and provide therapeutic targets for the cooperative treatment of BC.

To further decode the synergistic mechanism between Radix Sophorae Flavescentis and Rhizoma Heterosmilacis, we investigated the change in common target proportion under the condition of different degrees as thresholds. The results showed that following the increase in degree, the proportion of common target retention maintained a high growth trend. The proportion of common targets increased from 32.31 to 100% when the degree threshold increased from 1 to 28. Additionally, when the degree threshold increased from 1 to 28, the proportion of unique targets of Rhizoma Heterosmilacis decreased from 26.76 to 0%; however, when the degree threshold increased from 1 to 8, the proportion of specific targets of Radix Sophorae Flavescentis decreased from 41.02 to 0%. The results indicated that the common targets of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis have relatively high degrees. A higher degree indicates that these targets have a higher influence among all targets, and it also suggests that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis exert major synergistic effects by targeting common targets, which have a greater influence on the intervention response network.

To further decode the synergistic mechanism between them, we calculated the dynamic ratio change between common targets and unique targets of Radix Sophorae Flavescentis and Rhizoma Heterosmilacis using degree as the threshold. As shown in Figure 3, the proportion of common targets showed an overall upward trend from 32.31% (degree = 1) to 100% (degree = 28). However, the proportion of unique targets of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis revealed a gradual downward trend. The ratio of unique targets of Radix Sophorae Flavescentis changed from 41.02 to 0% with a degree ranging from 1 to 28, while the ratio of unique targets of Rhizoma Heterosmilacis changed from 26.67 to 0% with a degree range from 1 to 7.

FIGURE 3

Dynamic ratio changes between common targets and unique targets of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis.

In the C-T network, nodes with a higher degree usually represent the importance of all nodes. The results indicated that too much background noise may obscure the synergistic mechanism. Therefore, a model should be used to extract the core components.

Synergy Contribution Degree Calculation and Effect Verification

It is worth noting that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis may play their roles through common targets based on the C-T network. However, the core component groups and mechanisms of synergy require further elucidation. To solve these problems, an SCD calculation method was designed that considers the synergy of both the topology of the components in the C-T network and the dose of the components. The SCD values of the active components of CKI are shown in Supplementary Table S6.

According to the calculation results, the top two components with an SCD sum of 49.05% were oxymatrine (KS80) and matrine (KS20). Twenty-four components contribute to the effects of CKI on BC, with a total of 95.03%. Surprisingly, 24 components accounted for only 21.23% of all active components in CKI and could cover 631 targets (80.89% of all CKI targets) (Figure 4A). Thus, we define these 24 components as the high contribution group of synergistic components (HCGSC). These results show that HCGSC plays key roles in all active components of CKI and may clarify why the herbs in CKI generate synergistic and combination effects on BC.

FIGURE 4

(A) Cumulative radar chart of the ratio between component accumulation targets and the targets of CKI. The label within the chart represented the ratio (%) between component accumulation targets and the targets of CKI. (B) The dynamic ratio changes between common targets of HCGSC and unique targets of Radix Sophorae Flavescentis, and Rhizoma Heterosmilacis of HCGSC. (C) The network with targets and their enriched pathways (the degrees of the targets were higher than the third quartile of the targets in the C-T network of HCGSC) and their enriched pathways.

To evaluate the reliability of HCGSC, we defined two references for further comparison: the first reference was 70 genes, which overlapped with DEGs and CKI targets, and the second reference was the 17 pathways, which are overlapped by enriched DEGs pathways and CKI targets. The targets of HCGSC contained 62 genes, which could cover 88.57% (Supplementary Figure S4A) of the first reference, while the HCGSC targets enriched 158 pathways, which could cover 94.11% (Supplementary Figure S4B) of the second reference. Of these 62 overlapped genes, 39 were targeted by Rhizoma Heterosmilacis, while 43 were targeted by Radix Sophorae Flavescentis. These results confirm the reliability and accuracy of the HCGSC selection model. Two networks were constructed to analyze HCGSC. The first is a C-T network based on 24 components and 631 targets. The second is the target–pathway network based on 532 targets and 158 enriched pathways.

In the first network, the degree of targets ranged from 1 to 16, with the third quartile equal to 3. The change in the proportion of common targets under different conditions showed that the common targets had a higher degree (Figure 4B). The proportion of common targets increased from 27.26% to 100%. Additionally, with the increase in the threshold of the degree, the proportion of unique targets of Rhizoma Heterosmilacis decreased from 32.17 to 0%, while the proportion of specific targets of Radix Sophorae Flavescentis decreased from 40.57 to 0%. In the second network, we assigned the degree from the first network to the targets; targets with a degree higher than 3 (the third quartile of the targets in the C-T network) and their enriched pathways were selected to construct the network (Figure 4C). Visualization showed that the size of the target was positively related to the degree. We noted that the common targets (in blue diamond) displayed more connections with pathways and larger sizes than the others. However, the unique targets of Rhizoma Heterosmilacis (yellow triangle) and Radix Sophorae Flavescentis (yellow arrow) not only showed fewer connections with pathways, but were also smaller than the common targets. These two networks indicated that the common targets of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis in HCGSC have relatively high degrees, which shows that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis have a synergistic effect on the nodes with higher influence.

Potential Synergistic Mechanism Analysis of CKI in BC Treatment Based on HCGSC

GO Enrichment Analysis for CKI Based on HCGSC

To further interpret the potential synergistic mechanisms of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis, we performed GO-BP enrichment analysis for the targets of Radix Sophorae Flavescentis and Rhizoma Heterosmilacis in HCGSC. Radix Sophorae Flavescentis targets in HCGSC were enriched in 1517 GO-BP terms, while Rhizoma Heterosmilacis targets in HCGSC were enriched in 1394 GO-BP terms. A total of 927 joint GO-BP terms were found between them (Supplementary Figure S4C).

We selected the GO-BP terms with an FDR-adjusted p value lower than that of the first quartile of all GO-BP terms to build a network (Figure 5A). We found that 86.49% of the selected GO-BP terms belong to the commonly enriched GO-BP terms between the targets of Radix Sophorae Flavescentis and Rhizoma Heterosmilacis in HCGSC. This suggests that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis have potential synergistic mechanisms in the treatment of BC through these highly reliable, commonly enriched GO-BP terms based on HCGSC.

FIGURE 5

Components were represented in green. The GO-BP terms and pathways were represented in yellow. The size of the pathways was related to the FDR adjust p value. The lower FDR adjust p value represented the larger size. (A) The network of selected GO-BP terms based on the first quartile value of all GO-BP terms FDR adjust p value. The left of the network was the commonly enriched GO-BP terms. The upper right of the network was the unique GO-BP terms of Rhizoma Heterosmilacis. The lower right of the network was the unique GO-BP terms of Radix Sophorae Flavescentis. (B) The network of pathways between Rhizoma Heterosmilacis and Radix Sophorae Flavescentis. The center of the network was the commonly enriched pathways. The right of the network was the unique pathways of Rhizoma Heterosmilacis. The left of the network was the unique pathways of Radix Sophorae Flavescentis.

It is worth noting that most of the joint GO-BP terms were closely related to BC, such as the steroid metabolic process (GO:0008202, FDR adjusted p value = 2.24E-27) and steroid biosynthetic process (GO:0006694, FDR adjusted p value = 1.23E-18). Regarding the genes of these GO terms, estrone (E1) can protect women from breast cancer, while estradiol (E2) and estriol (E3) may enhance the risk (Lemon et al., 1966; Cohn et al., 2017). The expression of ATGL is correlated with tumor aggressiveness in vivo (Wang et al., 2017), which is related to the lipid catabolic process (GO:0016042, FDR adjust p value = 1.46E-22), lipid transport (GO:0006869, FDR adjusted p value = 4.04E-17), and lipid localization (GO:0010876, FDR adjusted p value = 1.16E-16). Including TRPCs, TPRVs, TRPMs, TRPA1, TRPPs, and TRPMLs (Berridge et al., 2003; Rohacs, 2005; Monteith et al., 2007; Venkatachalam and Montell, 2007; Doherty et al., 2015; Vangeel and Voets, 2019), calcium homeostasis plays an important role in the occurrence, development, and metastasis of breast cancer, which is related to the response to metal ions (GO:0010038, FDR adjust p value = 9.17E-13).

Additionally, the other GO-BP terms of targets for Radix Sophorae Flavescentis and Rhizoma Heterosmilacis in HCGSC, respectively, are also related to BC treatment. The ERK/MAPK and JAK2/PI3K signaling cascades in breast cancer cells can be activated by α7-nAChR activation (Chen et al., 2006; Nishioka et al., 2011; Kalantari-Dehaghi et al., 2015), while α9-nAChR overexpression is observed in tumor tissues compared with adjacent normal tissues (Lee et al., 2010). These genes are involved in synaptic transmission, cholinergic (GO:0007271, FDR adjust p value = 3.54E-19), and acetylcholine receptor signaling pathways (GO:0095500, FDR adjust p value = 1.91E-10) of Radix Sophorae Flavescentis in HCGSC. For the remaining GO-BP terms of Rhizoma Heterosmilacis in HCGSC, some steroid hormones, such as vitamin D, may have anticancer properties, while others may favor cancer progression, including estrogens and androgens (Restrepo-Angulo et al., 2020), which are related to regulation of steroid metabolic processes (GO:0019218, FDR adjusted p value = 2.49E-13).

Pathway Analyses Exploring Therapeutic Mechanisms of CKI Based on HCGSC

To further dissect the potential synergistic mechanisms of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis, we performed the KEGG pathway enrichment analysis for the targets of Radix Sophorae Flavescentis and Rhizoma Heterosmilacis in HCGSC, respectively. Targets of Radix Sophorae Flavescentis in HCGSC were enriched in 55 pathways, while those of Rhizoma Heterosmilacis were enriched in 54 pathways. Thirty-four commonly enriched pathways were found between them (Supplementary Figure S4D).

We used pathways to build a network (Figure 5B). Surprisingly, 45.33% of pathways belonged to the commonly enriched pathways between targets of Radix Sophorae Flavescentis and Rhizoma Heterosmilacis in HCGSC. This suggests that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis have potential synergistic mechanisms in the treatment of BC through these commonly enriched pathways based on HCGSC.

As shown in Figure 5B, most of the commonly enriched pathways were closely related to BC. Arachidonic acid metabolism (hsa00590, FDR adjusted p value = 4.47E-10), as a metabolic process, and the arachidonic acid (AA) pathway play key roles in carcinogenesis (Yarla et al., 2016). Additionally, PLA2s, COXs, LOXs, CYP-dependent monooxygenases, and their metabolites are known to play key roles in carcinogenesis (Tong et al., 2002; Go et al., 2015; Wu et al., 2015; Yarla et al., 2016). The migration and invasion of MDA-MB-231 breast cancer cells may be induced by linoleic acid (Serna-Marquez et al., 2017), which is related to linoleic acid metabolism (hsa00591, FDR adjusted p value = 4.47E-10). The PI3K-Akt signaling pathway (hsa04151, FDR adjusted p value = 2.78E-07) plays an important role in the tumorigenesis of breast cancer (Yuan and Cantley, 2008), and its activation may promote tumor progression in mice (Mei et al., 2018).

Additionally, other pathways of targets for Radix Sophorae Flavescentis and Rhizoma Heterosmilacis in HCGSC, respectively, are also related to BC treatment. For other pathways of Radix Sophorae Flavescentis in HCGSC, cellular Ca2+ signals have been implicated in the induction of apoptosis and regulation of apoptotic pathways (Berridge et al., 1998; Sergeev, 2005), which is related to the calcium signaling pathway (hsa04020, FDR adjusted p value = 1.43E-04). Regarding other pathways of Rhizoma Heterosmilacis in HCGSC, the higher alpha-linolenic acid content may reduce the risk of breast cancer, which is related to alpha-linolenic acid metabolism (hsa00592, FDR adjusted p value = 3.13E-06).

To further explain the synergistic mechanisms, we obtained 34 commonly enriched pathways between targets for Radix Sophorae Flavescentis SC and Rhizoma Heterosmilacis in HCGSC and 25 DEG-enriched pathways. By analyzing these pathways, we found that five pathways overlapped between 34 and 25 pathways, which we considered highly reliable pathways. These include the PI3K-Akt signaling pathway (hsa04151), cAMP signaling pathway (hsa04024), neuroactive ligand–receptor interaction (hsa04080), dopaminergic synapse (hsa04728), and ABC transporters (hsa02010). The PI3K-Akt and cAMP signaling pathways are reportedly related to BC, with a greater number of published studies (Mei et al., 2018; Dong et al., 2015). Next, we constructed a comprehensive pathway (Figure 6), including the PI3K-Akt and cAMP signaling pathways, using CyKEGGParser to explore the synergetic mechanisms of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis in the treatment of BC using CKI (Nersisyan et al., 2014).

FIGURE 6

Distribution of target proteins of CKI on the compressed BC pathway. The colors of the blanks represent different types of targets.

In the PI3K-Akt signaling pathway, Radix Sophorae Flavescentis acts on targets of the upstream pathway, such as Ras, GBY, and Syk, while Rhizoma Heterosmilacis acts on downstream targets, such as PKCs, SGK, Bcl-2, and Bcl-xL. These results indicate that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis have synergistic and complementary effects on cell survival, cell cycle progression, metabolism, survival signal, growth and proliferation, actin reorganization, vesicle transport, and glucose uptake. Additionally, in the cAMP signaling pathway, Rhizoma Heterosmilacis acts as a target of the upstream pathway, including Gi and PKA, while Radix Sophorae Flavescentis acts as a target of the downstream pathway, including FOS, NFKB, and CFTR, which are associated with pathways of hyperexcitability, cell death or survival (hippocampal neuron), and increased testicular AMH output (prepubertal Sertoli cells). As the pathogenic factors of BC are related to metabolism, cell survival, proliferation, cell death, and cell cycle (Jia et al., 2018), the above results suggest that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis can exert a synergistic effect on BC at multiple pathways.

Molecular Docking Validation of HCGSC in the Treatment of BC

To evaluate the function of HCGSC in BC treatment, molecular docking was employed to simulate the interaction between the ligands and the protein. Twenty-four components from the HCGSC of CKI and 85 proteins coded by 30 genes were used for molecular docking, and 18,254 docking relationships were obtained from the docking results (Supplementary Table S7). The result with the lowest affinity value in each component-protein was selected as the best docking relationship (2040 docking relationships), of which the affinity value ranged from −11.3 kcal/mol to −0.5 kcal/mol. Among the best docking relationships, three exhibited affinity values lower than −11.0 kcal/mol, BTL3 and BTL4 had the lowest affinity values of −11.3 kcal/mol, binding with protein 3ovv coded by PKA and protein 2w73 coded by CaM, respectively. BTL3 could bind 2w73 coded by CaM with the affinity values of −11.2 kcal/mol thereafter (Figures 7A–C). Based on the literature reports, improved binding between component–protein interactions should have a lower value of affinity (Elhenawy et al., 2019), and the histogram revealed that most of the results were concentrated in the lower affinity value position (Figure 7D).

FIGURE 7

(A) Docking result visualization of BTL3-3ovv (coded by PKA). (B) Docking result visualization of BTL4-2w73 (coded by CaM). (C) Docking result visualization of BTL3-2w73 (coded by Cam). (D) The frequency histogram of the affinity value. (E) The bar plot to display the composition of CDI and SCD in each TRI.

To reduce the influence of background noise of docking relationships with high affinity values, we selected the docking relationships that were lower than the average affinity value (−6.743) from all the best docking relationships as a highly reliable docking relationship (HRDR). Twenty-two components and 82 proteins coded by 29 genes were identified in HRDR with 1,255 docking relationships. Among the 22 components, 15 components were from Radix Sophorae Flavescentis, including KS1, KS2, KS3, KS4, KS11, KS16, KS19, KS20, KS52, KS53, KS70, KS77, KS90, KS91, and KS92, and seven components were from Rhizoma Heterosmilacis, including BTL1, BTL2, BTL3, BTL4, BTL8, BTL15, and BTL19. The 15 components of Radix Sophorae Flavescentis could bind 82 proteins coded by 29 genes, and the seven components of Rhizoma Heterosmilacis could bind 81 proteins coded by 29 genes. These results suggest that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis exert a synergistic effect on BC at multiple pathways.

To evaluate the therapeutic potential of components in HRDR, we designed a mathematical model to calculate the TRI of each component in the HRDR from the HCGSC. To further decode the mechanism of CKI in treating BC comprehensively, we integrated the CDI and SCD to calculate the TRI, which not only considered synergistic effects, but also considered the docking relationship. Detailed information on the TRI is listed in Supplementary Table S8. As shown in Figure 7E, KS20 (TRI = 1.08), KS91 (TRI = 1.01), BTL3 (TRI = 0.81), and KS1 (TRI = 0.78) were the top four components with TRI scores higher than 0.75. Among these four components, matrine (KS20), with the highest TRI, was confirmed by a number of published studies related to BC therapy (Li et al., 2015). The following three components with higher TRIs were selected for in vitro experiments to validate the reliability of our proposed strategy.

Luteolin and Trifolirhizin Inhibited Proliferation and Metastasis of MCF-7 Cells In Vitro

Based on the TRI results, trifolirhizin (KS91), beta-sitosterol (BTL3), and luteolin (KS1) were selected for in vitro experiments. The in vitro experiments proved that two of the three core components of HCGSC showed significant inhibitory effects on breast cancer cell proliferation and migration. The effects of luteolin and trifolirhizin on the viability of MCF-7 cells were determined using the MTT assay. As shown in Figures 8A,B, the proliferation of MCF-7 cells was suppressed by luteolin and trifolirhizin in a dose-dependent manner. Subsequently, the effects of luteolin and trifolirhizin on the migration of MCF-7 cells were investigated. As presented in Figures 7C,D, luteolin (40 and 80 Μm) and trifolirhizin (20 and 40 Μm) dramatically inhibited the migration of MCF-7 cells. These results indicated that luteolin and trifolirhizin markedly inhibited the proliferation and migration of MCF-7 cells.

FIGURE 8

Inhibitory effects of luteolin (A) and trifolirhizin (B) on the proliferation of MCF-7 cells. Inhibitory effects of luteolin (C) and trifolirhizin (D) on migration in MCF-7 cells. Data are represented as mean ± SEM (n = 3). * p < 0.05, ** p < 0.01, *** p < 0.001 as compared with that of the control group.

Discussion

Breast cancer is one of the most common malignant tumors among women worldwide, with the second highest death rate in female cancer patients (Siegel et al., 2020). At present, research regarding cancer treatment using TCM is still in the exploratory stage, and CKI has a wide range of applications in the treatment of tumors. However, there are few studies on CKI in the treatment of BC through systematic pharmacology. Therefore, we designed a systematic pharmacology strategy combining differentially expressed gene analysis, pharmacokinetics synthesis screening, target identification, synergy contribution degree, network analysis, therapeutic response index calculation, and experimental validation and provided a reference for this new method.

In this study, we examined the synergistic effect of CKI on BC at three levels. At the first level, the target coincidence of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis in CKI indicates that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis may act synergistically to exert a therapeutic effect in treating BC. First, we constructed a C-T network and found that most targets with a higher degree were mapped on the joint targets between Rhizoma Heterosmilacis and Radix Sophorae Flavescentis. These results suggest that Rhizoma Heterosmilacis and Radix Sophorae Flavescentis may act synergistically to treat BC via common targets and provide therapeutic targets for the cooperative treatment of BC. However, in the complex C-T network, we know neither the key synergistic components in treating BC nor their synergistic mechanisms.

To solve this issue, at the second level, we designed an SCD model with the advantages of reflecting the topology character of the C-T network and the component dose, which implemented a semi-quantitative system pharmacology model. Through the analysis of the SCD model, 24 components compose HCGSC, which can contribute to the effects of CKI on BC with a sum of 95.03% and map 80.89% (631/780) of all targets of CKI. Two references were employed to confirm the reliability and accuracy of the HCGSC selection model, including genes and pathways. The targets of HCGSC contained 62 genes, which could cover 88.57% of the first reference, while the HCGSC targets enriched 81 pathways, which could cover 70.59% of the second reference. Third, we performed GO-BP and KEGG enrichment analysis for Radix Sophorae Flavescentis and Rhizoma Heterosmilacis targets of HCGSC. We found that most of the overlapped GO-BP terms and pathways between them affect BC treatment, and most of them have a lower FDR adjusted p value.

At the third level, to further understand the synergistic mechanisms of CKI and infer the function of synergistic components and verify their effectiveness, three aspects were employed. First, a comprehensive pathway was constructed to explore the synergetic mechanism of Rhizoma Heterosmilacis and Radix Sophorae Flavescentis in the treatment of BC using CKI. In the PI3K-Akt and cAMP signaling pathways, Rhizoma Heterosmilacis and Radix Sophorae Flavescentis can activate upstream and downstream targets to play their role in the treatment of BC. These results suggest that CKI can produce a combined effect on BC. Second, we designed a TRI model, which provides the results of SCD and docking, to comprehensively evaluate the synergy contribution rate of the components in HCGSC and the binding ability of the components and proteins. TRI calculation results showed that four components displayed a TRI score higher than 0.75, including matrine (KS20, TRI = 1.08), trifolirhizin (KS91, TRI = 1.01), beta-sitosterol (BTL3, TRI = 0.81), and luteolin (KS1, TRI = 0.78), which may be the core components in the treatment of BC with synergistic mechanisms. Third, MTT and Transwell assays were employed to determine the function of luteolin and trifolirhizin. The results showed that they markedly inhibited the proliferation and migration of MCF-7 cells.

In this study, we proposed two novel mathematical models for the speculation of synergistic components and mechanisms. The method provided a methodological reference for decoding the synergistic mechanism of TCM in treating complex diseases. However, several limitations must be mentioned. 1) The precise synergistic mechanisms based on the calculation predictions warrant further validation. 2) Because the interaction between herbs is extremely complex, it may be incomplete to explore the synergistic mechanisms based on our model. Multiple synergistic effects combining with other techniques should be further considered. 3) The TRI model is proposed based on network topology methodology, which is more suitable for the case that herbs have a complex network. This algorithm used undirected network, which ignores the activation or inhibition effects of the disease targets.