Network Pharmacology and Molecular Docking Analysis on Molecular Targets and Mechanisms of Bushen Hugu Decoction in the Treatment of Malignant Tumor Bone Metastases

Purpose To explore the active compounds of the Chinese medicine prescriptions of Bushen Hugu Decoction (BHD) and demonstrate its mechanisms against malignant tumor bone metastasis (BM) through network pharmacology and molecular docking analysis. Methods The main components and targets of BHD were retrieved from the TCMSP database, and the targets were normalized by UniProt. The Herbs-Components-Targets network of BHD was established by Cytoscape. The main BM targets were obtained from GeneCards, TTD, DrugBank, and OMIM. STRING and Cytoscape were used to construct a PPI network and obtain hub genes. DAVID and Metascape were used for GO and KEGG enrichment analyses. According to the network topology parameters, the top 4 components were selected for molecular docking verification with the core targets. Results Compound–target network of BHD mainly contained 51 compounds and 259 corresponding targets including 107 BHD-BM targets. PPI interaction network and subnetworks identified ten hub genes. GO enrichment analysis found 1970 terms (p < 0.05), and 164 signaling pathways (p < 0.05) were found in KEGG, including PI3K-Akt signaling pathway, proteoglycans in cancer, prostate cancer, MAPK signaling pathway, and IL-17 signaling pathway. Molecular docking analysis showed that the active components of BHD, quercetin, luteolin, kaempferol, and aureusidin have good binding activity to the core targets. Conclusion The potential molecular target and signaling pathways were found for BHD major active components. It provides guidance for the future mechanism research of the BHD in malignant tumor bone metastasis. This study also established the foundation for the new strategy for the pharmacology study of Chinese medicine.


Introduction
Cancer incidence increased worldwide. Cancer, especially malignant tumor, ranks as the disease with the highest fatality rate after cardiovascular disease globally [1]. Most malignant tumor patients do not die from the primary disease directly but from complications such as cardiovascular disease, infection, secondary primary tumor, and metastases [2]. In recent years, with the development of preventive medicine, early diagnosis, and therapeutic regimens, the survival rate of cancer patients has improved. From 1999 to 2015, the overall cancer death rates decreased by 1.8% per year among men and by 1.4% per year among women [3]. However, the treatment for patients at the advanced stage with multisite metastasis is still very challenging. Metastasis is a common cause of death in malignant tumors. Bone is one of the most common sites of malignant tumor metastatic after the lung and liver. Bone metastases occur in approximately 70% of prostate cancer or breast cancer patients and 30-40% of lung cancer patients [4][5][6]. Bone metastases can lead to skeletal-related events (SREs), including hypercalcemia, pathological fractures, spinal cord compression, and    [7][8][9][10]. The standard cares for bone metastases in cancer patients are bisphosphonates and denosumab. Bisphosphonates, the inhibitors of osteoclast activity, can inhibit tumor growth by blocking prenylation. However, clinical studies showed minimal efficiency on tumors due to limited uptake at tumoral sites [11]. Denosumab is a monoclonal antibody against the receptor activator of nuclear factor-kappa-Β ligand (RANKL). The RANK/RANKL/OPG axis is critical for the regulation of bone formation and resorption. Inhibition of RANKL may serve as an ideal therapeutic strategy for excessive bone resorption [12]. However, whether long-term use of denosumab has side effects on bone tissues and the immune system is still unknown, further research is needed [13]. Standard care combined with radiation has become a new strategy for the treatment of BM. Unfortunately, this treatment strategy is only partially effective and is palliative, not curative [14,15].
Chinese medicine is a traditional treatment method original to China. A variety of Chinese medicine prescriptions, herbs, and extracts have been confirmed with the potential to treat malignant tumors by experimental studies and clinical trials [16]. Combining ancient wisdom and modern medical science, Chinese medicine has become an essential part of the comprehensive treatment of malignant tumors in China and several other Asian countries [17] [18,19].
Cistanches Herba extract phenylethanol glycosides increase the sensitivity of liver cancer patients to oxaliplatin by inhibiting HIF-1α signaling pathway [20].
Cuscutae Semen, Dipsaci Radix, Drynariae Rhizoma, and Eucommiae Cortex can inhibit bone destruction through the immune system and inflammation-related pathways [21][22][23][24]. In this study, we used network pharmacology technology to investigate the potential mechanism of BHD on bone metastases and molecular docking technology to screen the major effective components and targets for BHD.  [25]. TCMSP is built on the framework of systems pharmacology for herbal medicines to ensure the confidence of results. The toxic pharmacokinetic property parameters were set as follows: oral availability ðOBÞ ≥ 30%, drug-likeness ðDLÞ ≥ 0:18, and half-life ðHLÞ ≥ 4 [26]. The filtered active components were used to screen their protein targets. The targets of active compounds were searched through PubMed, CNKI, and Wanfang databases. Protein target gene names were canonicalized through the UniProt database (https://www.uniprot .org). "Herbs-components-targets" network of BHD was constructed based on network topology parameters (NTP) through Cytoscape 3.7.2. The larger the node, the more targets associated with it, and the more likely it will be a key ingredient in BHD.

Construction of the BHD-BM Gene PPI Interaction
Network. The bone metastases-related targets were identified in the GeneCards database (https://www.genecards.org), OMIM database (http://www.omim.org), TTD database (http://bidd.nus.edu.sg/group/cjttd), and DrugBank database (https://www.Drugbank.ca) using "Bone Metastasis", "Bone Metastases", "Osseous Metastasis", "Bone Metastasis of Malignant Tumor" and "Metastatic Carcinoma of Bone" as keywords. The targets were supplemented through literature search. Use Venn script in R language to get the inter-section of the targets of BHD components and BM. BHD-BM genes were uploaded to STRING11.0 (https://String-db .org) to build a PPI interaction network. Initial screening of the BHD-BM gene by setting the degree score. The data were corrected by setting the biological species to 'Homo sapiens'. The PPI interaction network of BHD-BM genes was visualized with Cytoscape using 0.9 as the highest confidence level. Topological analysis of BHD-BM genes by Cytohubba in Cytoscape to obtain the hub genes [27].

GO and KEGG Enrichment
Analyses. GO and KEGG enrichment analyses of BHD-BM genes were performed using DAVID (https://david.ncifcrf.gov/) online analysis tool and Metascape platform (http://metascape.org/gp/ index.html). A p < 0:05 was used to screen BHD-BMrelated biological process (BP), molecular function (MF), and cellular component (CC). Then, the top ten results were chosen for visualization. Similarly, a p < 0:05 was used to  5 BioMed Research International screen KEGG pathway analysis. BHD compounds, BHD-BM genes, and top 20 model pathways were used in KEGG enrichment analysis. The NTP (degree, betweenness, and closeness) was analyzed by the built-in tool of Cytoscape to build the BHD compound, BM target, and pathway network.
2.4. Active Component-Target Docking. Molecular docking is a well-established computer-based structural method for structurally docking small molecules to target proteins and assessing their binding affinity [28]. It is widely used in new drug discovery. Negative docking energy indicates that the small molecule can bind effectively and autonomously to the target protein [29]. The lower the binding energy, the tighter the bond, and the more effective it is. In this study, the core components of BHD and genes were subjected to bulk molecular docking. The 3D structures (SDF format) of potentially active compounds were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) [30]. And save it as PDBQT format through AutoDock Tools 1.5.6. The 3D structures (PDB format) of the targets were obtained from the RCSB-PDB database (http://www .RCSB.org/) [31]. Use AutoDock Tools to remove water molecules, separate proteins, add nonpolar hydrogen, calculate the Gasteiger charge of the structure, and save it as a PDBQT file. The targets were used as receptors and set to be rigid. The compounds were used as ligands and set to be flexible. AutoDock Vina 1.1.2 was used to perform bulk molecular. The active sites for docking were determined by adjusting the original ligands on the different receptors by setting up a grid box size of 15 × 15 × 15 points with a spacing of 1.0 Å between grid points, covering almost the entire range of favourable protein binding sites. Each set of ligands and receptors produced 9 conformations. Collect the binding energy of the best-affinity conformation for a heatmap. The pocket structures and specific docking conformations of top 3 ligands bound to top 3 targets were visualized by PyMol 2.5 [32]. Table 1, sixty-two active ingredients of BHD were found in TCMSP database, including Cuscutae Semen (CS, 11 compounds), Rehmanniae Radix Praeparata (RRP, 2 compounds), Eucommiae Cortex (EC, 16 compounds), Dipsaci Radix (DRD, 6 compounds), Cistanches Herba (CH, 5 compounds), Drynariae Rhizoma (DR, 16 compounds), Sinapis Semen (SS, 3 compounds), and Radix Angelicae Biseratae (RAB, 3 compounds). A total of 51 active ingredients were obtained after weight removal (Table 1). A total of 752 targets were obtained for the active compounds of BHD: CS (344 targets), RRP (34 targets), EC (156 targets), DRD (18 targets), CH (22 targets), DR (141 targets), SS (20 targets), and RAB (17 targets). After removing duplications, 259 targets were obtained ( Figure 1). Quercetin, luteolin, kaempferol, and aureusidin were identified as the key components of BHD.

GO and KEGG Enrichment
Analyses. GO enrichment analysis found 1,970 terms. Multiple targets of BHD were closely related to BM. The major biological process enriched for BHD-BM has a response to an inorganic substance, positive regulation of cell death, response to cytokine, positive regulation of cell migration, and regulation of apoptotic signaling pathway. The molecular function related to BHD-BM has kinase binding, protein kinase binding, ubiquitin-like protein ligase binding, DNA-binding transcription factor binding, and cytokine receptor binding. BHD-BM-related cellular components have membrane raft, transcription regulator complex, vesicle lumen, serine/threonine protein kinase complex, and extracellular matrix (Figure 3(a)). KEGG pathway enrichment found 164 terms. The main pathway relative to HBD-BM included PI3K-Akt signaling pathway, proteoglycans in cancer, MAPK signaling pathway, IL-17 signaling pathway, and microRNAs in cancer ( Table 2). Figure 3(b) shows the BHD compound, BM target, and pathway network constructed by Cytoscape using the top 20 signaling pathways.

Component-Target Docking Analysis.
According to the NTP analysis, the core compounds (Table 3) (Table 4) PTGS2,  HSP90AA1, RELA, MAPK1, AKT1, MAPK3, JUN,  MAPK14, CHUK, PRKACA, TNF, NFKBIA, TP53, IL6, CASP3, and VEGFA. The molecular docking results showed that the binding energies of the core components to the core genes were both negative, indicating that they can bind autonomously and tightly (Figure 4(a). Visualize the molecular docking results of the top 3 core compounds with the top 3 core genes by PyMol (Figures 4(b)-4(d)). Of these, quercetin was the most closely related to PTGS2, with five hydrogen bonds formed with the PTGS2 amino acid residues GLY-536, ASN-375, LEU-145, TYR-373, and ASN-375.

Discussion
Bone is one of the common sites of cancer metastasis. The unique structure and the dynamic metabolism pattern in bone make bone metastasis a different mode of malignant tumor metastasis. The mechanism of BM is still not very clear yet.
The "seed and soil" theory for tumor metastasis proposed by Paget in 1889 believed that the spread of tumor cells is determined by the interaction between cancer cells (seed) and host organs (soil) [33]. After more than a century's research, it is well recognized that the occurrence, development, and metastasis of tumors are not linear but involve multiple parallel overlapping routes [34,35]. Bone  hsa04115  p53 signaling pathway  17 -25.120  BAX, CCND1, BCL2, BCL2L1, CASP3, CASP8, CASP9, CDK1, CDK2, CDK4,  CDKN1A, IGFBP3, MDM2, SERPINE1, PTEN, TP53, and CHEK2    11 BioMed Research International microenvironment (BME) and tumor microenvironment (TME) are important factors for bone metastasis [36]. BME is composed by bone osteoblast cell, osteoclast cell, osteogenic cell, osteocyte cell, chondrocytes, adipocytes, fibroblasts, hematopoietic and mesenchymal stem cells, bone matrix, and various growth factors IGFs, BMPs, TGFβ1, PDGFs, etc. [37]. TME includes immune cells, blood vessels, extracellular matrix (ECM), fibroblasts, lymphocytes, bone marrow-derived inflammatory cells, and signaling molecules. Interactions between tumor cells and nontumor cells in TME play important roles in cancer development, progression, and metastasis [38]. Bone marrow is an important hematopoietic organ and immune organ. In addition, immune cells have multiple interactions with bone marrow stromal cells [39]. The crosstalk between BME and TME may promote the metastatic potential of tumor cells.
The interaction between TME and BME involves multiple pathways as a comprehensive network. However, during the drug development process, it is very common to identify or design pharmacologically effective agents specifically targeting a single molecule or site. Drugs that act on a single molecular target often show unsatisfactory therapeutic effects or have large side effects when used to treat complex diseases in clinical. Therefore, drug development has gradually shifted from a "single target and single drug" model to a "network target multicomponent therapy" model [40,41]. Multicomponents contained in Chinese medicine prescriptions (CMP) together with their targets constitute a huge drug interaction network. Due to this reason, CMP maybe not as good as the chemical drugs in the efficiency, which brings challenges to the evidence-based medicine verification of CMP. Network pharmacology is a bioinformatic network construction and network topology analysis strategy based on high-throughput data analysis, virtual computing, and network database retrieval [42]. It provides the possibility to study the complex network of CMP components-targets.
In this study, a total of 51 active components and 259 targets for BHD and 1086 BM targets were obtained. Combination of these two dataset, 107 BHD-BM intersection targets, was obtained finally. PPI interaction network analysis obtained the hub targets. GO and KEGG enrichment analyses showed that multiple targets of active components of BHD were closely related to BM. Membrane raft, transcription regulator complex, vesicle lumen, serine/threonine protein kinase complex, and extracellular matrix are the top cellular components found in GO and KEGG enrichment analyses. These cellular components are closely related to the formation of extracellular vesicles (EVs). EVs include exosomes, activation-or apoptosis-induced macrovesicles, and apoptotic bodies, which are key players in intercellular communication [43]. EVs are a heterogeneous group of membrane-limited vesicles loaded with lipids, various proteins, and multiple genetic materials such as DNAs, mRNAs, and noncoding RNAs [44]. Myeloid-derived suppressor cells (MDSC) are a heterogeneous population of immature myeloid cells with immunosuppressive activity. Tumor cells and stromal cells regulate the activation and expansion of MDSCs in the TME through EVs, thereby affecting the body's antitumor immune response [45].
KEGG enrichment analysis of BHD-BM showed that BHD compounds can affect BM through various pathways including regulation of apoptosis, immunity, inflammatory response, and treatment of primary tumors. A3 (quercetin), A2 (kaempferol), DR13 (aureusidin), DR9 (luteolin), A1 (beta-sitosterol), EC13 (syringetin), DR10 (xanthogalenol), and B2 (stigmasterol), the top 8 active components of BHD, were obtained after the bioinformatic analysis. Their combination to core targets was verified by molecular docking, which confirmed the multitarget and multichannel therapeutic effect of BHD on BM. PI3K-Akt signaling pathway was identified as one of the top signaling pathways in BHD-BM network. PI3K-Akt could participate in multiple cellular events to promote tumor initiation, progression, and metastasis. PI3K/Akt significantly increases bone metastases and osteolytic bone lesions in prostate cancer by mediating the stabilization of histone methyltransferase WHSC1 [46]. TLRs are key pathogensensing receptors to promote development and activation of the immune system. In TME, TLRs are widely expressed in the tumor-infiltrating immune cells and have critical roles in the prometastasis inflammatory microenvironment [47]. The activated TLR2 : TLR6 can induce TNF-alpha secretion by myeloid cells and enhance lung cancer bone metastasis growth [48]. Notably, IL-17 signaling pathway, T cell receptor signaling pathway, Th17 cell differentiation, and osteoclast differentiation enriched for BHD-BM were closely related to bone marrow-mediated immune activity. Homeostasis of osteoblasts and osteoclasts maintains bone stability. IL-17 is a marker for Th17 cells, a special subset of CD4 + helper T cells [49]. Th17 cells are positively associated with osteoclast formation in patients with myeloma [50]. By blocking the differentiation and function of Th17 cells, bone marrow mesenchymal stem cells can improve multiple sclerosis in an experimental autoimmune encephalomyelitis model [51]. All this evidence suggested that BHD may regulate bone marrow-mediated immunity by affecting Th17 cells.
In this study, we built the network between BHD compounds and BM target and screened and verified the in vivo combination of BHD components and their target employing network pharmacology and molecular docking analysis. Network pharmacology showed that quercetin and kaempferol were the top 2 core compounds of BHD.
Quercetin and kaempferol, as natural bioflavonoids, are widely found in a large number of natural herbs and have broad-spectrum antitumor effects [52,53].

Conclusion
This study established the foundation of the new strategy for the pharmacology study of Chinese medicine. This study 12 BioMed Research International shows that the main active components of BHD can participate in the regulation of tumor and bone microenvironment through extracellular vesicles and affect the occurrence and development of malignant tumor bone metastases through various mechanisms such as immunity, apoptosis, and inflammation. It provides a basis for the study of the molecular mechanism of BHD and provides a new option for the treatment of bone metastases.

Data Availability
All the data can be obtained from the open-source website we provide, and the conclusion can be drawn through the analysis of the relevant software.

Conflicts of Interest
The authors declare no competing interests.