Network Pharmacology Prediction: The Possible Mechanisms of Cinobufotalin against Osteosarcoma.

OBJECTIVE: To explore the active compounds and targets of cinobufotalin (huachansu) compared with the osteosarcoma genes to obtain the potential therapeutic targets and pharmacological mechanisms of action of cinobufotalin on osteosarcoma through network pharmacology.
METHODS: The composition of cinobufotalin was searched by literature retrieval, and the target was selected from the CTD and TCMSP databases. The osteosarcoma genes, found from the GeneCards, OMIM, and other databases, were compared with the cinobufotalin targets to obtain potential therapeutic targets. The protein-protein interaction (PPI) network of potential therapeutic targets, constructed through the STRING database, was inputted into Cytoscape software to calculate the hub genes, using the NetworkAnalyzer. The hub genes were inputted into the Kaplan-Meier Plotter online database for exploring the survival curve. Functional enrichment analysis was identified using the DAVID database.
RESULTS: 28 main active compounds of cinobufotalin were explored, including bufalin, adenosine, oleic acid, and cinobufagin. 128 potential therapeutic targets on osteosarcoma are confirmed among 184 therapeutic targets form cinobufotalin. The hub genes included TP53, ACTB, AKT1, MYC, CASP3, JUN, TNF, VEGFA, HSP90AA1, and STAT3. Among the hub genes, TP53, ACTB, MYC, TNF, VEGFA, and STAT3 affect the patient survival prognosis of sarcoma. Through function enrichment analysis, it is found that the main mechanisms of cinobufotalin on osteosarcoma include promoting sarcoma apoptosis, regulating the cell cycle, and inhibiting proliferation and differentiation.
CONCLUSION: The possible mechanisms of cinobufotalin against osteosarcoma are preliminarily predicted through network pharmacology, and further experiments are needed to prove these predictions.

1. Introduction

Osteosarcoma (OS) [1] is the most common cancer in orthopedics, which occurs in adolescents or over-65-year-old patients. The most common lesions of osteosarcoma are the actively growing metaphysis, including the distal femur and proximal tibia, which can seriously affect the patient in daily life and activity. Pain and swelling [2] are the most common symptoms of osteosarcoma in the early stage. The pain is mostly intermittent at first and is often confused with growing pain, leading to a later diagnosis. Surgery [3] is the main treatment for osteosarcoma at present. Adjuvant chemotherapy and neoadjuvant chemotherapy can significantly improve the prognosis of patients with nonmetastatic osteosarcoma. With the current comprehensive treatment, close to 2/3 of patients suffering from osteosarcoma can be cured; also, the limb rate reaches 90%-95%.

Traditional Chinese medicine (TCM) has been treating tumors and cancers for more than two thousand years. Simiao Sun, a famous physician in the Tang Dynasty, supported that the skin of toads eliminates malignant swelling. Cinobufotalin [4], extracted from the skin of the Chinese giant salamander or black sable through modern technology, has significant detoxification, swelling, and pain relief effects; is used for advanced tumors; and has become a class A reimbursed antitumor drug by China Medical Insurance. At present, there are few studies on the antiosteosarcoma of cinobufotalin. Therefore, we used network pharmacology to analyze the possible targets and pathways of cinobufotalin against osteosarcoma to point out the direction of the research on the mechanism of cinobufotalin and provide more theoretical and data support for the modernization of TCM.

2. Materials and Methods

2.1. The Active Compound and Target Predictions of Cinobufotalin

Taking “huachansu,” “dry scalp,” “cinobufotalin,” and “cinobufacin” as keywords, we searched the possible compounds of cinobufotalin through PubMed (https://pubmed.ncbi.nlm.nih.gov/), CNKI (https://www.cnki.net/), Wanfang (http://www.wanfangdata.com.cn/index.html), VIP (http://www.cqvip.com/), and other literature databases. The literature search screened 90 compounds. Through the PubChem database (https://http://pubchem.ncbi.nlm.nih.gov/), we obtained 28 active ingredient structures. We explored the targets of cinobufotalin via the Comparative Toxicogenomics Database (CTD, http://ctd.mdibl.org/) and Traditional Chinese Medicine Database and Analysis Platform (TCMSP, https://tcmsp-e.com/). Active compounds and therapeutic targets were retrieved and inputted to Cytoscape software (version 3.7.2) to construct the “active-compound-target” network of cinobufotalin.

2.2. Osteosarcoma Gene Prediction and Potential Therapeutic Target Screening

With “osteosarcoma” as the keyword, we explored the GeneCards database (https://www.genecards.org/), DisGeNET database (https://www.disgenet.org/), Online Mendelian Inheritance in Man (OMIM, https://omim.org/), and Therapeutic Target Database (TTD, http://db.idrblab.net/ttd/) to obtain osteosarcoma genes, as well as combined the genes of those databases as osteosarcoma genes. The targets of cinobufotalin and genes of osteosarcoma were imported into the Venn production website [5] (http://www.bioinformatics.com.cn/) to obtain the potential therapeutic targets between cinobufotalin and osteosarcoma.

2.3. PPI Interaction Network and Survival Prognosis of Hub Genes

The common targets were imported into the STRING database (https://string-db.org/) to construct a protein-protein interaction (PPI) network with a confidence level ≥ 0.400. And then, the PPI file was imported into Cytoscape 3.7.2 to construct a PPI network diagram. NetworkAnalyzer was used to calculate the nodes' degree, and the top 10 targets were regarded with degree as the hub genes, which included TP53, ACTB, AKT1, MYC, CASP3, JUN, TNF, VEGFA, HSP90AA1, and STAT3. The relationship between the mRNA expression levels of hub genes and the prognosis (survival rate) was analyzed through the Kaplan-Meier Plotter online database (http://kmplot.com/analysis/).

2.4. Functional Enrichment Analysis

Under the condition of P value < 0.05, the potential therapeutic targets were enriched and analyzed through the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/), including the Cellular Component (CC), Molecular Function (MF), Biological Process (BP) of Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The possible mechanisms of action of cinobufotalin against osteosarcoma were analyzed to extract the 5 smallest P value terms of three types in GO and 15 KEGG pathways for analysis.

3. Results

Through literature retrieval and the PubChem database, 28 active compounds of cinobufotalin were selected. For the active compounds, CTD and TCMSP databases were used to search for therapeutic targets. By importing the compounds into the CTD and TCMSP databases, we found that there are relatively many targets in the CTD database. For instance, uracil, adenosine, thymidine, adenine, and other components which exist in all organisms have hundreds of targets in the CTD database. However, only a few targets can be explored in the database for bufalin, cinobufagin, bufotenine, and other unique ingredients in cinobufotalin. This situation can lead to unreality of predicted targets and pathways. In contrast, there are few targets in the TCMSP database. As for the unique ingredients of cinobufotalin, the TCMSP database hardly includes their targets, which may also lead to untrue prediction data. After discussion by the research team, according to the distribution of the obtained active ingredients in nature, we set the TCMSP database of the widely existing ingredients in the active ingredients as the target prediction database and the CTD database of the unique ingredients in cinobufotalin. From this, a total of 184 targets were selected, and the “active-compound-target” network diagram constructed by Cytoscape 3.7.2 is shown in Figure 1. The network contains 28 active compounds and 184 targets, with a total of 265 edges. On average, each active compound regulates 9.46 targets. The active compounds with the most regulated targets are bufalin, adenosine, oleic acid, cinobufagin, and gamabufotalin, with degrees of 60, 43, 23, 21, and 19, respectively.

Figure 1

“Active-compound-target” map of cinobufotalin. Yellow nodes represent the active compounds; green nodes represent therapeutic targets.

3544 osteosarcoma genes were confirmed through the GeneCards database, 2283 genes from the DisGeNET database, 11 genes through the OMIM database, and 7 genes from the TTD database. After deleting the duplicates, 4188 targets were retained, and the cinobufotalin targets and osteosarcoma genes were imported into the Venn map production website, in which 128 potential therapeutic targets were obtained. The Venn map is shown in Figure 2.

Figure 2

Venn diagram of potential therapeutic targets of cinobufotalin against osteosarcoma. Green circle represents cinobufotalin and blue circle represents osteosarcoma. Their intersection represents potential therapeutic targets.

We imported the potential therapeutic targets into the STRING database to construct the PPI network. Since 2 genes are not directly related to other targets, the network contains a total of 126 targets. And then, we exported the data to construct the network in Cytoscape 3.7.2. Calculating the degree through NetworkAnalyzer, we intuitively and clearly revised the node and edge font size for the network diagram (Figure 3). The targets with the largest degree are TP53 (degree = 103), ACTB (degree = 100), AKT1 (degree = 99), MYC (degree = 90), CAPS3 (degree = 88), and JUN (degree = 88).

Figure 3

PPI network of cinobufotalin in the treatment of osteosarcoma. Green nodes represent potential therapeutic targets of cinobufotalin against osteosarcoma.

The hub genes by Cytoscape 3.7.2 were introduced into the Kaplan-Meier Plotter online database, and the survival curve of sarcoma was analyzed. Correlation analysis showed that TP53, ACTB, MYC, TNF, VEGFA, and STAT3 were correlated with the median survival time of sarcoma (P < 0.05). The median survival time of the TP53, MYC, and VEGFA low expression group is better than that of the high expression group, as shown in Figures 4(a), 4(c), and 4(e); the median survival time of ACTB, TNF, and STAT3 high expression group is better than that of the low expression group, as shown in Figures 4(b), 4(d), and 4(f).

Figure 4

Survival analysis of hub genes: (a) for TP53, (b) for ACTB, (c) for MYC, (d) for TNF, (e) for VEGFA, and (f) for STAT3. Red curve represents high expression group, and black curve represents low expression group.

Function enrichment analysis obtains 666 GO terms, including 515 Biological Process terms, which mainly include negative regulation of apoptotic process, response to drug, apoptotic process, activation of cysteine-type endopeptidase activity involved in apoptotic process, and response to hypoxia. 67 Cellular Component terms include cytosol, nucleus, nucleoplasm, death-inducing signaling complex, cytoplasm, and cytosol, and 84 terms of Molecular Function mainly include identical protein binding, protein binding, enzyme binding, protein kinase binding, and transcription factor binding. 118 KEGG pathway terms include pathways in cancer, Hepatitis B, apoptosis, proteoglycans in cancer, TNF signaling pathway, colorectal cancer, toxoplasmosis, Chagas disease (American trypanosomiasis), small cell lung cancer, and NF-kappa B signaling pathway. The details are shown in Tables 1 and 2.

Table 1

GO analysis of the potential therapeutic targets.

Serial number	Term	Count	P value	FDR
GO:0043066	Negative regulation of apoptotic process	32	2.31E-21	4.37E-18
GO:0042493	Response to drug	25	4.37E-18	4.12E-15
GO:0006915	Apoptotic process	30	1.22E-16	7.70E-14
GO:0006919	Activation of cysteine-type endopeptidase activity involved in apoptotic process	15	1.42E-15	6.70E-13
GO:0001666	Response to hypoxia	17	1.95E-13	7.37E-11
GO:0005829	Cytosol	70	2.11E-20	4.95E-18
GO:0005634	Nucleus	76	3.27E-12	3.82E-10
GO:0005654	Nucleoplasm	51	2.48E-11	1.93E-09
GO:0031264	Death-inducing signaling complex	6	3.03E-10	1.77E-08
GO:0005737	Cytoplasm	67	1.71E-08	7.76E-07
GO:0042802	Identical protein binding	37	7.55E-20	2.74E-17
GO:0005515	Protein binding	113	6.19E-19	1.12E-16
GO:0019899	Enzyme binding	26	2.41E-18	2.91E-16
GO:0019901	Protein kinase binding	22	5.87E-13	5.32E-11
GO:0008134	Transcription factor binding	19	3.73E-12	2.71E-10

GO: Gene Ontology; FDR: false discovery rate.

Table 2

KEGG pathway of the potential therapeutic targets.

Serial number	Term	Count	P value	FDR
hsa05200	Pathways in cancer	52	8.80E-34	8.18E-32
hsa05161	Hepatitis B	35	7.47E-31	3.47E-29
hsa04210	Apoptosis	25	1.04E-27	3.23E-26
hsa05205	Proteoglycans in cancer	29	6.79E-19	1.58E-17
hsa04668	TNF signaling pathway	23	1.18E-18	2.20E-17
hsa05210	Colorectal cancer	19	2.58E-18	4.00E-17
hsa05145	Toxoplasmosis	22	3.92E-17	5.20E-16
hsa05142	Chagas disease (American trypanosomiasis)	21	2.03E-16	2.36E-15
hsa05222	Small cell lung cancer	19	1.32E-15	1.36E-14
hsa04064	NF-kappa B signaling pathway	18	3.49E-14	3.25E-13
hsa04620	Toll-like receptor signaling pathway	19	7.66E-14	6.42E-13
hsa05152	Tuberculosis	23	8.29E-14	6.42E-13
hsa04932	Nonalcoholic fatty liver disease (NAFLD)	21	3.75E-13	2.69E-12
hsa05133	Pertussis	16	8.68E-13	5.77E-12
hsa05219	Bladder cancer	13	1.40E-12	8.68E-12
hsa04115	p53 signaling pathway	15	2.78E-12	1.62E-11
hsa04621	NOD-like receptor signaling pathway	14	4.14E-12	2.25E-11
hsa05160	Hepatitis C	19	4.36E-12	2.25E-11
hsa05164	Influenza A	21	5.62E-12	2.75E-11
hsa04722	Neurotrophin signaling pathway	18	8.52E-12	3.96E-11

KEGG: Kyoto Encyclopedia of Genes and Genomes; FDR: false discovery rate.

4. Discussion

Network pharmacology [6, 7] is based on systems biology, integrating pharmacology, high-throughput sequencing, genomics, and other technologies, emphasizing multicomponent, multipathway, and multitarget regulation mode. It is an important means to predict the mechanism of action of herbs and TCM preparation, providing theoretical basis and future directions for exploring the treatment of diseases with TCM. This study explored the potential targets of the main active ingredients of cinobufotalin through network pharmacology and explored the possible mechanism of cinobufotalin against osteosarcoma.

Through literature and database predictions, we have determined that the main active compounds of cinobufotalin include bufalin, adenosine, oleic acid, cinobufagin, and gamabufotalin. Bufalin, cinobufagin, and gamabufotalin are unique ingredients of cinobufotalin. Bufalin [8, 9] is a cardiotonic steroid exerting antitumor activity in a variety of cancers by inducing apoptosis and anti-inflammatory pathways. In vivo research [10] had found that bufalin can alleviate the pain caused by sarcoma in bone models. In bone destruction, research [11, 12] showed that bufalin promotes osteosarcoma cell apoptosis by regulating miR-221, endogenous, and exogenous apoptotic pathways. Cinobufagin [13] increase the suppressed transient outward potassium currents of dorsal ganglion cells in a rat model of cancer pain and bone pain to exert analgesic effects. In vitro [14-16], cinobufagin can pass a variety of pathways, including the exogenous pathway, Notch signaling pathway, and MAPK signaling pathway to promote osteosarcoma cell apoptosis. Compared with the effect of cisplatin alone, the combination of low-dose cisplatin significantly inhibits cell viability and motility, induces apoptosis and cell cycle arrest in the S phase, and inhibits tumor growth and metastasis, as well as prolongs the exposure of mice. Beside, longer survival time [17] in the OS xenograft model significantly inhibited the Notch signaling pathway. Gamabufotalin [18] inhibits the viability of osteosarcoma cells and tumorigenesis by blocking the TGF-β/periostin/PI3K/AKT signaling pathway.

By analyzing the PPI network of the potential therapeutic targets, hub genes such as TP53, ACTB, AKT1, MYC, CASP3, JUN, and TNF were selected. TP53 is the most frequently changed tumor suppressor gene in osteosarcoma. Mutations [19] in the P53 reduce the 2-year survival rate of osteosarcoma, and the targeted mutation of TP53 may be a potential strategy for the treatment of osteosarcoma [20]. In more than 10% of osteosarcoma patients [21], MYC is overexpressed and promotes cell invasion by activating MEK-ERK pathways. In vitro inhibition of MYC [22] can inhibit the migration, clone formation, and spheroid growth of osteosarcoma cells. ACTB, the promoter of β-actin, as a component of the cytoskeleton, is an important mediator of intracellular movement. AKT1 is a key protein in the PI3K/AKT signaling pathway. Lots of researches [18, 23] have shown that bufalin and gamabufotalin can promote osteosarcoma apoptosis through the PI3K/AKT signaling pathway. CASP3 is the most important terminal shear in the process of apoptosis. Cinobufosin [16, 24, 25] can induce the apoptosis of osteosarcoma cells by activating caspase-3. JUN is a major protein kinase involved in the composition of the MAPK signaling pathway, playing an important role in the proliferation, differentiation, and apoptosis of osteoblasts. JUN inhibitors [26] can reduce the cell proliferation and metastasis of osteosarcoma. Epidemiological studies [27] confirmed that the proinflammatory factor TNF-α is associated with an increased risk of osteosarcoma, and TNF can regulate the survival of osteosarcoma cell through a variety of ways [28, 29]. In this study, bufalin regulates almost all the hub genes, cinobufagin regulates MYC and CASP3, and gamabufotalin regulates MYC, AKT1, and VEGFA. By analyzing the median survival time of hub genes to sarcoma, we found that the downregulation of TP53, MYC, and VEGFA and the upregulation of ACTB, TNF, and STAT3 can increase the median survival time of patients. But there are still few research reports, and further experimental proof is needed.

Through enrichment analysis, it is explored that the main way of cinobufotalin in the treatment of osteosarcoma is to promote cell apoptosis and inhibit tumor cell proliferation. A large number of studies [14, 15, 18, 30] have found that cinobufotalin can pass through mitochondrial apoptosis, Notch signaling pathway, MAPK signaling pathway, NF-kappa B signaling pathway, and PI3K/AKT signaling pathways to treat osteosarcoma.

5. Conclusion

Through network pharmacology, we explored the possible mechanism of cinobufotalin against osteosarcoma, using bioinformatics methods to obtain 28 active ingredients and 184 potential targets, of which 128 may be related to osteosarcoma. TP53, MYC, VEGFA, ACTB, TNF, and STAT3 are considered to be hub genes, and adjusting these targets can increase the median survival time of patients. Function enrichment analysis showed that the possible mechanism of cinobufotalin in the treatment of osteosarcoma is to promote sarcoma cell apoptosis, regulate tumor cell cycle, and inhibit cell proliferation and differentiation. At present, research on the pathogenesis and treatment of osteosarcoma is focused on apoptosis and autophagy [31]. Therefore, the network pharmacology of cinobufotalin against osteosarcoma is consistent with the current research understanding, but further experiments are needed to prove these predictions.