Screening of target genes of CFF-1 active ingredients in prostate cancer
The flowchart of this study was shown in Fig. 1. According to the screening results of the TCMSP database, CFF-1 contained 112 active ingredients, and it's found that licorice had the most target genes (Table S1). In table S1, the same target genes could be regulated by different small ingredients from different Chinese medicine, hence the regulation of CFF-1 on prostate cancer was diverse. According to five databases (GeneCards, OMIM, PharmGkb, TTD, and DrugBank), the prostate cancer target genes were obtained (Fig. S1a). After crossing the target genes of the above analysis methods, 206 target genes regulated by CFF-1 were obtained in prostate cancer (Fig. S1b and Table S2).
Biofunctions and pathways of target genes of CFF-1 active ingredients in prostate cancer
In this study, GO and KEGG enrichment analysis was employed to figure out the functions and pathways of the target genes of CFF-1 active ingredients in prostate cancer. According to the GO enrichment analysis of 206 genes, it was found that functions of these target genes in response to metal ion, membrane raft, and nuclear receptor activity were the most significant in BP (Biological Process), CC (Cellular Component), and MF (Molecular Function), respectively (Fig. 2a). In the results of KEGG enrichment analysis, prostate cancer signaling pathways including IL-17 signaling pathway, EGFR tyrosine kinase inhibitor resistance, and PI3K-Akt signaling pathway, etc. were obtained (Fig. 2b). Besides, other cancer-related pathways, such as colorectal cancer, bladder cancer, small cell lung cancer, hepatocellular carcinoma, and chronic myeloid leukemia (Fig. 2b), were also significantly enriched. Therefore, these target genes (marked in red) played important roles in cancer pathways (including estrogen signaling pathway and MAPK signaling pathway) (Fig. 2c).
Construction of the merge network of core PPI and drug active ingredients
Based on the screening results of the active ingredient and 206 target genes in the TCMSP database, we constructed a network of active ingredients and target genes (Fig. 3). In this network, many target genes, such as PTGS2, ESR1, MAPK14, AR, etc. could be regulated by a variety of active ingredients of CFF-1. According to the STRING database, the PPI network of 206 target genes was obtained (Fig. 4a). The core PPI network was obtained from the PPI network of 206 target genes by analyzing the relationship between nodes in the PPI network of 206 target genes (Fig. 4b). In the core PPI network, the genes, such as TP53, AKT1, MAPK family genes, ESR1, etc. had complex regulatory relationships, indicating that inhibiting these core genes might affect the entire PPI network in prostate cancer. After conjoint analysis of the core PPI network and the network of active ingredients and target genes, the merge network was rebuilt (Fig. 5). In the merge network, ESR1 and MAPK14 were repeatedly regulated by CFF-1 active ingredients, indicating that these two genes played major roles in the process of CFF-1 inhibiting prostate cancer.
Potential ligands targeting ESR1 and MAPK14 were obtained by virtual screening
According to the results of Do_virtual_screening (virtual screening tool developed by us, https://github.com/daizao/Do_virtual_screening), we found that Glabrene and (2S)-2-[4-hydroxy-3-(3-methylbut-2-enyl)phenyl]-8-8-dimethyl-2-3-dihydropyrano[2-3-f]chromen-4-one had the smallest binding affinities for ESR1 and MAPK14 among 46 ligands (Table S3). The docking result of Glabrene and ESR1 was displayed in 2D and 3D, hydrogen bond interaction between Arg and ligand was found (Fig. 6a). Docking result of 2S)-2-[4-hydroxy-3-(3-methylbut-2-enyl)phenyl]-8-8-dimethyl-2-3-dihydropyrano[2-3-f]chromen-4-one and MAPK14 was displayed in 2D and 3D, the hydrophobic interaction played a major role (Fig. 6b).
Molecular dynamics simulation of the binding of ESR1 and MAPK14 to ligands
Based on the results of the two ligands and protein complexes, molecular dynamics simulations were performed. The RMSD (root-mean-square deviation) of the carbon α skeleton represents the stability of the three-dimensional structure model of the protein-ligand complex during molecular dynamics simulation. After 100ns molecular dynamics simulation of ESR1 and glycyrrhizin, the RMSD value balanced at 0.31nm (Fig. 7a). The short-range Coulomb force and the short-range Lanna-Jones potential were − 85 (kJ/mol) and − 182 (kJ/mol) (Fig. 8a) at 100ns, respectively. Hence the complex of ESR1 and ligand was dominated by electrostatic force and Van der Waals force. After 150ns molecular dynamics simulation of MPAK14 and (2S)-2-[4-hydroxy-3-(3-methylbut-2-enyl)phenyl]-8-8-dimethyl-2-3-dihydropyrano[2-3-f]chromen-4-one, the RMSD value balanced at 0.53nm (Fig. 7b). The short-range Coulomb force and the short-range Lanna-Jones potential were − 6 (kJ/mol) and − 223 (kJ/mol) (Fig. 8b) at 150ns, respectively. Hence the complex of MAPK14 and ligand was dominated by Van der Waals force.