Network pharmacology and experimental evaluation strategies to decipher the underlying pharmacological mechanism of Traditional Chinese Medicine CFF-1 against prostate cancer

Active compounds in CFF-1 and target screening of CFF-1 on PCa

Initially, 78 active compounds in CFF-1, each with OB ≥30% and DL ≥0.1 in CFF-1, were collected from TCMSP platform. Subsequently, 70 active compounds with p < 0.05, score >20 were obtained from BATMAN-TCM platform. After integrating these findings, 112 active compounds were selected for further study (Table 1). The targets for the candidate compounds in CFF-1 were explored from TCMSP and BATMAN-TCM, identifying 131 and 848 putative targets, respectively. There were 41 overlapping targets among the two sets. Ultimately, 938 targets for CFF-1 active components were acquired by integrating the overlapping targets (Figure 1A). Additionally, 2682 target genes related to PCa were acquired from the Genecards by setting the correlation score >20 and 495 target genes were acquired from OMIM databases. A total of 3022 target genes acquired from the two sets after eliminating the overlaps (Figure 1B). Nine hundred and thirty-eight targets for drug active components and 3022 targets for PCa were screened out by Perl language program and R language software. In total, 359 overlapping target genes were recognized as key targets related to both CFF-1 and PCa for further analyses (Figure 1C).

Figure 1. Target screening of CFF-1 on PCa. (A) Venn diagrams showing CFF-1 targets obtained from TCMSP and BATMAN-TCM. (B) Targets related to PCa acquired from Genecards and OMIM. (C) The intersection of targets for both CFF-1 and PCa.

Compound-target network and analysis

Upon entering the above 359 key targets into STRING, we obtained a key targets PPI network of CFF-1 on PCa. Subsequent analysis of this PPI network focused on “degree” was applied to select the target in the core position (Figure 2A). The top 30 target genes with high degree were shown in Figure 2B. The 53 core targets were further screened by setting interaction score ≥0.9 and degree ≥20. Cytoscape software was applied to construct the herbs-active ingredients-target network, containing 112 nodes (55 for candidate active ingredients and 53 for core targets) and 253 edges. In the network, the nodes with more edges might play essential roles in the pharmacological processes, and 12 nodes (edge ≥5) were determined as the main targets of CFF-1 on PCa for further analysis, including NCOA2, RXRA, ESR1, NCOA1, PPARG, IL1B, TNF, IKBKB, NR3C1, IL4, IL6 and PRKCA (Figure 2C).

Figure 2. Compound-target network and analysis. (A) Key targets PPI network of CFF-1 on PCa. (B) Bar graph showing the top 30 targets with high degree. (C) The herbs-active ingredients-target network of CFF-1 on PCa.

Molecular docking verification

As previously mentioned, PIK3R1, AKT1, MAPK1, MAPK3 and SRC were the top five hub targets in PPI network. These targets were then docked with their respective active ingredients. The ligand-receptor binding energy values, indicative of binding stability, were presented in Figure 3A. Generally, a binding energy lower than −5 kcal/mol is considered indicative of stable binding. Notably, coryneine establishes two hydrogen bonds with GLU-17 and THR-18 in PIK3R1, while baicalein forms three hydrogen bonds with GLU-91, HIS-89 and HIS-13 in AKT1. Diosgenin forms one hydrogen bond with GLY-16 in AKT1. Coryneine and digitalis glycoside bound to MAPK1 with binding energy values of −5.6 kcal/mol and −9 kcal/mol, respectively. Among these, the binding affinity of digitalis glycoside to MAPK3 was the strongest, with a value of −12.9 kcal/mol. Additionally, nerolidol forms one hydrogen bond with GLU40 in SRC (Figure 3B–3H).

Figure 3. Representative images of molecular docking. (A) The results of ligand-receptor binding energy values. (B) Coryneine-PIK3R1. (C) Baicalein-AKT1. (D) Diosgenin-AKT1. (E) Coryneine-MAPK1. (F) Digitalis glycoside-MAPK1. (G) Digitalis glycoside-MAPK3. (H) Nerolidol-SRC.

Enrichment analysis

To ascertain the involved pathways of CFF-1 on PCa, we carried out KEGG pathway analysis of 359 key targets. A total of 172 involved pathways were identified (Supplementary Table 1). Figure 4A displayed the top 20 enriched pathways, suggesting that the anti-cancer effect of CFF-1 on PCa likely results from a relatively complex and multi-pathway synergistic effect.

Figure 4. KEGG pathway enrichment and GO biological process analysis of key targets of CFF-1 on PCa. (A) The top 20 enriched pathways of KEGG pathway enrichment analysis. (B) The top 20 enriched terms of GO biological process analysis.

To further explore the biological role of involved targets of CFF-1 against PCa, the GO process analysis of 359 key targets was performed. A total of 238 GO terms were found (Supplementary Table 2). Figure 4B displayed the top 20 enriched terms related to PCa, suggesting that CFF-1 may exhibit its therapeutic effects through the above biological processes.

Prognostic model construction based on key targets

A total of 359 key target genes were subjected to LASSO regression analysis to identify genes for construction of risk score model for PCa patients in the TCGA database. The risk score for predicting DFS of each PCa patient was calculated as follows: Risk score = NFATC1 × 0.234 + ARG1 × 0.001 + DNMT3B × 0.315 − NR3C1 × 0.025 + RXRB × 0.105 + HPRT1 × 0.033 + SI × 0.037. A 7-gene signature was constructed. Patients were classified into high- and low-risk groups based on the median risk score (Figure 5A). A risk curve and a scatter plot were applied to show the risk score and the survival status of each PCa patient, respectively. Most recurrent cases were distributed in the high-risk group (Figure 5B). The expression profile of candidate genes indicated that NFATC1, ARG1, DNMT3B, NR3C1, RXRB, HPRT1 and SI were highly expressed in the high-risk group, except for NR3C1 (Figure 5C). PCa patients in the high-risk group had a significantly worse DFS than those in the low-risk group. AUC of ROC curve for the DFS prediction of risk score model was 0.867 (Figure 5D, 5E).

Figure 5. Risk score for target gene signature and outcome in PCa patients. (A) Risk score of a 7-gene signature for predicting DFS. (B) Disease status and duration of cases. (C) Heatmap of the 7 gene expression in PCa patients. (D) Kaplan-Meier curve for DFS in the low- and high-risk groups. (E) ROC curve for the DFS prediction of risk score model. (F) Risk score of a 11-gene signature for predicting PFS. (G) Progression status and duration of cases. (H) Heatmap of the 11 gene expression in PCa patients. (I) Kaplan-Meier curve for PFS in the low- and high-risk groups. (J) ROC curve for the PFS prediction of risk score model.

The risk score for predicting PFS was calculated as follows: Risk score =ALDOA × 0.0003 + DNMT3B × 0.222 + CSF2 × 0.099 + DNMT1 × 0.031 + EZH2 × 0.054 + ARRB2 × 0.019 − ESRRG × 0.025 + APOA2 × 0.007 + RXRB × 0.076 + HPRT1 × 0.021 − RDH11 × 0.0001. A 11-gene signature was constructed (Figure 5F). Progression patients were mainly distributed in the high-risk group (Figure 5G). The expression profile showed that ALDOA, DNMT3B, CSF2, DNMT1, EZH2, ARRB2, APOA2, RXRB and HPRT1 were highly expressed in the high-risk group, except for ESRRG and RDH11 (Figure 5H). Kaplan–Meier curves showed that PFS was significantly worse in high-risk patients than low-risk patients. AUC of ROC curve for the PFS prediction was 0.773 (Figure 5I, 5J).

For OS, though few patients reached the endpoint (10 of 495 patients), the risk score model was also constructed. The risk score for predicting OS was calculated as follows: Risk score = AR × 0.048 + AGTR1 × 0.062 + PPARD × 0.427 + PHB × 0.090 + RPS6KB1 × 0.620 + FADD × 0.821 − DNMT1 × 0.667 +AKR1C3 × 0.075. An 8-gene signature was constructed (Supplementary Figure 1A). More dead cases were found in the high-risk group (Supplementary Figure 1B). The heat map suggested that AR, AGTR1, PPARD, PHB, RPS6KB1, FADD, DNMT1 and AKR1C3 were overexpressed in the high-risk group, except for DNMT1 (Supplementary Figure 1C). Patients in the high-risk group had significantly shorter OS compared with those in the low-risk group. AUC of ROC curve for the OS prediction was 0.992 (Supplementary Figure 1D, 1E). Taken together, these gene signatures might effectively predict the prognosis of PCa and may act as potential targets for PCa therapy.

We then explored possible associations of risk scores with ESTIMATE score using the ESTIMATE algorithm. Using the CIBERSORT algorithm, the correlation between risk score and the infiltration of 22 immune cell subtypes were assessed. The risk score for predicting OS of PCa patient was negatively associated with stromal score, immune score, and ESTIMATE score, while positively with the infiltration of mast cells resting (Supplementary Figure 2A). The risk score for predicting DFS was negatively related to the infiltration of plasma cells, while positively Tregs and macrophages M2 (Supplementary Figure 2B). The risk score for predicting PFS was positively related to stromal score, immune score, ESTIMATE score, the infiltration of Tregs macrophages M1 and macrophages M2, while negatively with plasma cells and mast cells resting (Supplementary Figure 2C).

DNMT3B, RXRB and HPRT1 were the common target genes affecting both PFS and DFS in PCa patients. GO functional annotations and KEGG pathways enrichment were applied to evaluate the biological significance of DNMT3B, RXRB and HPRT1 in PCa. The results indicated that DNMT3B, RXRB and HPRT1 were widely involved in immune regulation (Supplementary Figure 3).

CFF-1 inhibited PCa cells proliferation and induced apoptosis

To evaluate the effect of CFF-1 on PCa as postulated from network pharmacology analysis, clonogenic assay, CCK8 and Edu assays were conducted in PC-3 cell line treated with different concentrations of CFF-1 (0, 2, 4, 6, 8 and 10 mg/ml). The cell colony formation efficiency of PC-3 cells was decreased dose dependently compared to the negative control (Figure 6A, 6C). Furthermore, a dose-dependent reduction in proliferation active cells was observed after 24 hours of CFF-1 treatment in the EdU assay (Figure 6B, 6D). The CCK8 assay confirmed a dose-dependent decrease in PCa cell viability (Figure 6E). In addition, CFF-1 treatment increased the percentage of apoptotic PC-3 cells dose dependently (Figure 7A–7H).

Figure 6. CFF-1 inhibited PCa cells proliferation. (A) The colony formation efficiency of PC-3 cell lines with varying concentrations of CFF-1 (mg/mL). (B) The EdU assay of proliferation active cells after 24 hours of CFF-1 treatment. (C) Quantification of the colony formation efficiency with bar graph. (D) Quantitative results of the EdU assay with bar graph. (E) The cell viability was determined by CCK8 assay after CFF-1 administration. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001.

Figure 7. Apoptosis analysis of CFF-1-treated PC-3 cells by flow cytometry after 24 hours of treatment. Each panel corresponds to a different treatment condition: (A) NC, (B) DMSO, and (C–G) increasing concentrations of CFF-1 at 2, 4, 6, 8, and 10 mg/ml, respectively. (H) Quantitative analysis of the percentage of apoptotic cells across different treatment conditions. ^**p < 0.01, ^***p < 0.001.

CFF-1 inhibited PCa cells migration and invasion

Subsequently, wound healing and transwell experiments were conducted to explore the migration and invasion abilities of PC-3 cells treated with CFF-1. We found that fewer cells migrated to the scratch site after 24 hours of CFF-1 treatment dose dependently (Figure 8A, 8C). Furthermore, the migration and invasion ability of PC-3 cells were significantly inhibited after treatment with increasing CFF-1 (Figure 8B, 8D, 8E).

Figure 8. CFF-1 inhibited PCa cells migration and invasion. (A) Wound healing abilities of PC-3 cells in the presence of NC, DMSO and varying concentrations of CFF-1, with images captured at 0 and 24 hours post-wounding. (B) Representative images showing transwell migration and invasion assay of PC-3 cells treated with NC, DMSO and varying concentrations of CFF-1. (C) Quantitative analysis of the wound healing assay with different treatment conditions. (D) Quantitative analysis of the transwell migration assay. (E) Quantitative analysis of the transwell invasion assay. ^**p < 0.01, ^***p < 0.001.

CFF-1 attenuated proliferation pathways in PCa cells

According to network pharmacology analysis, the PI3K-Akt, HIF-1, TNF, EGFR-TKI resistance and PD-1 checkpoint pathway might play a crucial role in regulating PCa cell proliferation and survival by CFF-1. Then, we evaluated the expressions level of the common key targets of the five pathways. Pretreatment of CWR22Rv1 and PC3 cells with CFF-1 (2, 6 and 10 mg/ml) resulted in apparent repression of P-ERK1, NFκB1, RELA, P-mTOR, VEGFA, PD-L1, P-PI3K, P-AKT, TNF-α, P-EGFR and HIF-1α in dose dependently (Figure 9A, 9B). ELISA assays indicated that the secretory IL-6 of CWR22Rv1 and PC3 cells was blocked by CFF-1 dose dependently (Figure 9C, 9D). These findings suggested that the above five pathways might be crucial for CFF-1 anti-cancer effect in PCa.

Figure 9. The relative expressions of related proteins with CFF-1 treatment on PCa cells. (A) The relative expressions of related proteins after 24 hours of CFF-1 treatment on PC-3 cells. (B) The relative expressions of related proteins after 24 hours of CFF-1 treatment on 22RV1 cells. (C) The secretory levels IL-6 after 48 hours of CFF-1 treatment on PC-3 cells. (D) The secretory levels IL-6 after 48 hours of CFF-1 treatment on 22RV1 cells.

Screening for active ingredients of CFF-1

All ingredients of four main compounds (guizhi, fuzi, shudihuang, huangjing) in CFF-1 were acquired from TCM Systems Pharmacology Database and Analysis Platform (TCMSP, http://tcmspw.com/tcmsp.php) and a Bioinformation Analysis Tool for Molecular mechanism of Traditional Chinese Medicine (BATMAN-TCM, http://bionet.ncpsb.org/batman-tcm/). The criteria to screen for active ingredients were set as the oral bioavailability (OB) ≥30% with the drug similarity (DL) ≥0.1, and p < 0.05 with score >20, respectively.

Target prediction related to CFF-1 and PCa

The potential targets of active ingredients in CFF-1 were obtained from TCMSP and BATMAN-TCM. Targets related to PCa were acquired from Gene Cards database (https://www.genecards.org/) and Online Mendelian Inheritance in Man database (OMIM, http://www.omim.org/) with a keyword “prostate cancer”. The UniProt database (https://www.uniprot.org/) was applied to standardize all target information, and genes without UniProt ID from human samples were eliminated. Perl software was applied to screen intersecting key targets related to both CFF-1 and PCa.

Protein–protein interaction analysis and network construction

STRING database (https://string-db.org/) is used to construct the relationship between multiple proteins. We put the key targets of CFF-1 and PCa into the STRING database Version 11.0, selecting the Homo sapiens for the species to perform the PPI network with Cytoscape software Version 3.9.1. Targets with interaction score ≥0.9 and degree ≥20 were regarded as core targets. In order to study the relationship among CFF-1, the core targets and PCa, main herbs, active ingredients and core targets were imported into the Cytoscape software to construct the herbs-active ingredients-target network.

Molecular docking

Molecular docking was carried out to evaluate the interaction between the top five targets, which exhibited high degree, and their respective active ingredients. Using the PubChem database, we acquired the 3D molecule structures of ingredients. The protein crystal structures of the targets were retrieved from the RCSB PDB database. The PyMOL software facilitated the removal of water molecules and the separation of ligands. For the conversion of small molecules and target proteins into pdbqt format, AutoDock Tools 1.5.6 was utilized. Molecular docking was performed using AutoDock Vina software. Visualization and analysis of the docking results were achieved using the PyMOL tool.

Enrichment analyses

The Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment were carried out through gene set enrichment analysis (GSEA) in R Version 3.7.0. Biological significance was defined as p < 0.05, and Q < 0.05.

Prognostic model construction based on key targets

RNA sequencing datasets and survival information for PCa were obtained from TCGA data portal. LASSO Cox regression was performed to screen genes with prognostic values from the key targets using R package “glmnet”. A risk score formula for predicting prognosis was established as follows: risk score = expgene1 × βgene1 + expgene2 × βgene2 + … expgenen × βgene. Kaplan–Meier method was used to assess the prognostic significance of the risk score model on PCa by the R package “survival”. Receiver operating characteristic (ROC) curve analysis was performed to explore the predictive power of risk score using the R package “timeROC”.

We applied the CIBERSORT algorithm to assess the correlation between risk score and the infiltration of 22 immune cell subtypes [41]. Based on the ESTIMATE algorithm, we evaluated the relationship between risk score and immune score, stromal score and ESTIMATE score [42].

Cell culture and drug preparation

The human PCa cell lines CWR22Rv1 and PC-3 were cultured in RPMI-1640 medium, enriched with 10% fetal calf serum (Gibco, Waltham, MA, USA) and maintained in a humidified atmosphere with 5% CO₂ at 37°C. CFF-1, a TCM herbal mixture, has been detailed in our previously published study [11]. A final concentration of 0, 2, 4, 6, 8 and 10 mg/mL of CFF-1 was used to treat cells.

Cell viability assay

PCa cell viability was detected by Cell Counting Kit-8 (CCK8) assay kit (Beyotime, Guangzhou, China) based on the manufacturer’s instructions. PC-3 cells were seeded into 96-well plates at a density of 2 × 10³ cells/well. Each well was measured at 450 nm for its absorbance.

EdU assay

In order to further evaluate cell proliferation ability, PC-3 cells were measured using the BeyoClick™ EdU-555 assay according to the manufacturer’s protocol (Beyotime, Guangzhou, China). The signals were measured by fluorescence microscopy Olympus CKX53.

Clonogenic assay

PC-3 cells were transplanted into 6-well plates with a density of 2 × 10³ cells/well and then the cells were treated with or without CFF-1 incubated for 10 days at 37°C. Cells were fixed in 4% paraformaldehyde and stained with crystal violet. The numbers of visible colonies were counted under an optical microscope.

Apoptosis assay

PC-3 cells were treated with various concentrations of CFF-1. After 24 h, apoptosis was detected by the annexin V-FITC apoptosis detection kit (cat. #640932, Biolegend, USA) adhering to the manufacturer’s protocol.

Wound healing assay

PC-3 cells were added into 6-well plates and incubated until 100% confluence. After scraping in a straight line using a 200 μl pipette tip, the cells were rinsed thrice with PBS. Varying concentration of CFF-1 was administered to each well. The cell migration data were acquired with an inverted microscope Olympus IX51 at 0 and 24 h incubation and assessed using Image-Pro Plus 7.0 software.

Transwell migration and invasion assay

The transwell migration and invasion assay was performed using a transwell chamber system with an 8-μm pore polycarbonate membrane (Thermo Fisher Scientific, Waltham, MA, USA). For migration assay, 1 × 10⁵ cells containing different concentrations of CFF-1 intervention and medium supplemented with 2% serum were plated onto 24-well chambers. For invasion assay, we diluted the Matrigel 1:4 with serum-free medium and seeded it to the upper chambers. Then, the cells were placed onto the upper chambers. Both assays were conducted as previously described [43].

Western blot

Total protein was extracted using RIPA buffer containing proteinase inhibitor (Best Biological, Jiangsu, China). Western blot was conducted as previously described [12]. The primary antibodies for P-ERK1, ERK1, NF-κB1, RELA, P-mTOR, mTOR, VEGFA, PD-L1, P-P13K, P13K, P-AKT, AKT, TNF-α, P-EGFR, EGFR and HIF-1α were purchased from Bioss (Beijing, China). Total protein level was normalized to β-actin.

ELISA

PC3 and CWR22Rv1 cells treated with different concentration of CFF-1 were cultured for 48 h. The supernatant of cell culture was collected and centrifuged at 1,000 g for 20 min. ELISA kits (Mlbio, Shanghai, China) were used to examine levels of IL-6 following the manufacture’s instruction. We measured the absorbance at 450 nm with a microplate reader.

Statistical analysis

Statistical analyses of the in vitro PCa cell assays were conducted using SPSS 20.0 software. One-way univariate analysis of variance (ANOVA) was employed to analyze data obtained from at least three independent experiments, which are presented as the mean ± standard deviation. Statistical significance was defined as p < 0.05.

Availability of data and materials

All data included in this study are available by contacting the corresponding authors.