Network pharmacology reveals the potential mechanism of Baiying Qinghou decoction in treating laryngeal squamous cell carcinoma

Context: Baiying Qinghou as a traditional Chinese medicine decoction shows anticancer property on laryngeal squamous cell carcinoma. However, little is known about the precise mechanism of Baiying Qinghou detection against laryngeal squamous cell carcinoma. Objective: This study was aimed to explore potential mechanism of therapeutic actions of Baiying Qinghou decoction on laryngeal squamous cell carcinoma. Materials and Methods: The active chemical components of Baiying Qinghou decoction were predicted, followed by integrated analysis of network pharmacology and molecular docking approach. The network pharmacology approach included target protein prediction, protein-protein interaction network construction and functional enrichment analysis. Results: Sitosterol and quercetin were predicted to be the overlapped active ingredients among three Chinese herbs of Baiying Qinghou decoction. The target proteins were closely associated with response to chemical, response to drug related biological process and cancer related pathways such as PI3K-Akt signaling, HIF-1 signaling and Estrogen signaling pathway. The target proteins of TP53, EGFR, PTGS2, NOS3 and IL1B as the key nodes in PPI network were cross-validated, among which EGFR, IL1B, NOS3 and TP53 were significantly correlated with the prognosis of patients with laryngeal squamous cell carcinoma. Finally, the binding modes of EGFR, IL1B, NOS3 and TP53 with quercetin were visualized. Discussion and Conclusion: Quercetin of Baiying Qinghou decoction showed therapeutic effect against laryngeal squamous cell carcinoma by regulating TP53, EGFR, NOS3 and IL1B involved with drug resistance and PI3K-AKT signaling pathway. TP53, EGFR, NOS3 and IL1B may be the candidate targets for the treatment of laryngeal squamous cell carcinoma.


INTRODUCTION
Head and neck cancer (HNC) is the seventh most common malignant tumor worldwide and occurs in nasopharynx, larynx, and thyroid [1,2]. Head and neck squamous cell carcinoma (HNSCC) accounts for over 90% cases of head and neck cancer (HNC) and mainly derives from the mucosal surfaces of the upper aerodigestive tract [3]. Laryngeal squamous cell carcinoma (LSCC) represents the second most common histological subtypes of HNC with an increasing incidence rate [4]. Patients with LSCC are usually diagnosed at a late clinical stage and the survival rate is lower due to regional or distant metastases. Thus, LSCC has caused noticeable medical and economic burden worldwide.
Currently, several mainstream options against LSCC include surgery, radiotherapy, chemotherapy and chemo-radiotherapy [5]. However, the overall survival of patients with LSCC has not been remarkably improved because of chemotherapy or drug resistance and undesired effects. Traditional Chinese medicines (TCM) has been considered as an attractive alternative therapy for conquering cancers in China [6]. TCM focuses on restoring body balance and boosting immunity by the synergistic effects of various active ingredients [7]. Some Chinese herbs have been frequently utilized for the treatment of different diseases, such as cardiovascular diseases, diabetes and caners [8][9][10]. A recent research highlights that erchen plus huiyanzhuyu decoction can inhibit the growth of laryngeal carcinoma by modulating STAT3/cyclin D1 signaling pathway [11]. TCM exerts promising therapeutic effect on various malignant tumors. Network pharmacology is proposed by Hopkins et al. and aims to explore the multilevel interactions of diseases, genes, and drugs as a whole [12]. This systems pharmacology is based on systems biology, computational biology and omics theory to evaluate the therapeutic effects of Chinese medicines on several diseases [13]. Gao et al. indicate that 8 herbs regulate multiple hepatocellular carcinoma-related genes and are strongly associated with prognosis by a network pharmacology approach [14]. Huang et al. suggest Huanglian Jiedu decoction has therapeutic roles in cancers such as hepatocellular carcinoma by building the herb-compound, compound-protein, proteinpathway, and gene-disease networks [15]. Notably, Bai Ying Qing Hou (BYQH) decoction is a Chinese medicinal formula and widely used for the treatment of LSCC in China, which is primarily composed of five key traditional Chinese herbs: Solanum lyratum thumb (30 g), Scutellaria barbata (24 g), Duchesnea indica (24 g), Solanum nigrum (30 g), and Actinidia chinensis planch (30 g). However, the pharmacological mechanism of BYQH decoction in the treatment of LSCC has not been clarified.
In this study, we attempted to explore the mechanism of BYQH decoction in the treatment of LSCC by network pharmacology analysis. First, the key active chemical constituents and their target proteins of BYQH decoction were screened. Then, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of gene targets were carried out. Besides, protein-protein interaction (PPI) and target proteins were cross-validated followed by survival analysis. Finally, molecular docking was used to identify the binding mode of key active ingredients and target proteins of BYQH decoction.

Screening active constituents of BYQH decoction
The traditional Chinese medicine systems pharmacology database (TCMSP) is an herbal repository and provides the system information about Chinese herbal medicines [16]. The absorption, distribution, metabolism, and excretion (ADME) system has been successfully used to evaluate the pharmacokinetics characteristics of chemical compounds [17,18]. Herein, the chemical composition of five key components of BYQH decoction was firstly searched from TCMSP, mainly including the number of the pharmaceutical ingredient, molecule name, molecular weight, lipid-water partition coefficient, the number of hydrogen-bonding donor or acceptor, oral bioavailability (OB), blood brain barrier, Caco-2 permeability (Caco-2), drug-likeness (DL), and half-life (HL). OB ≥ 30% and Caco-2 ≥ −0.4 indicate good OB and permeability of molecules. The mean DL value of drugs collected in database is 0.18. Then, active chemical compounds were screened in ADME system according to the criteria of OB ≥ O30%, DL ≥ 0.18 and Caco-2 ≥ −0.4 according the previous description [19].

Screening protein targets of active chemical compounds
TCMSP and DrugBank (https://www.drugbank.ca/) databases can predict the relationships between drugs and corresponding targets. Moreover, DrugBank database provides the detailed information on experimental and investigational drugs [20]. The active chemical compounds were subjected to TCMSP and DrugBank databases to search their potential targets. Then, the target proteins were mapped to corresponding gene symbols using string database. For those proteins that did not map to any gene symbols, they were AGING manually retrieved by Universal Protein (UniProt) database to map to corresponding gene symbols. These potential gene targets were used for following analysis.

Functional enrichment analysis of target genes
To investigate the biological function of the target genes, the gene ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were carried out. There are three categories for GO terms, including biological process (BP), cellular component (CC) and molecular function (MF). Gene symbols were input into string database. The parameters were set as the default values. The significantly enriched GO and pathway terms with p < 0.05 were collected. Then, the top 10 GO terms in each category and top 10 significant pathways were visualized in bubble chart by using R language.

Target protein-protein interaction (PPI) analysis
The underlying relationships among target proteins were analyzed based on the string database according to the default parameters [21]. The Cytoscape software was used to construct and visualize the PPI network. Following this, the topological characteristics of the PPI network were also evaluated. The node degree (connectivity) was calculated.

Cross-validation of target proteins by the database retrieval
DisGeNET is a versatile platform that comprehensively collects human gene-disease associations (over 380,000 associations between more than 16,000 genes and 13,000 diseases) for the validation of computationally predicted genes for human diseases [22,23]. In this study, the genes related to LSCC were firstly searched from this database using "carcinoma of larynx" as the keyword. Phenopedia database also provides systematical genetic association studies and summarizes the information about the association between genes studied and a particular disease. Therefore, the known genes associated with LSCC were obtained by retrieving Phenopedia database with the keyword of "laryngeal neoplasms". Subsequently, the overlapping genes among target genes of BYQH decoction, LSCCrelated genes in DisGeNET or Phenopedia database were extracted by jvenn tool [24,25]. These genes may be associated with the therapeutic effects of BYQH decoction against LSCC.

Survival analysis of candidate genes
Firstly, expression profiles of HNSCC samples deposited in the Cancer Genome Atlas (TCGA) dataset were downloaded from UCSC Xena (https://xenabrowser.net/datapages/). Based on the clinical information, the HNSCC samples were selected according to following criteria: 1) samples were primary tumor tissues and the primary site was at the larynx; 2) the clinical information of samples with regard to TNM stage, clinical stage, grade and followup were complete; 3) the samples without survival time were deleted. Finally, a total of 105 laryngeal cancer samples were retained for the next prognostic analysis.
To investigate whether key candidate genes are associated with prognosis of patients with laryngeal cancer, the survival analyses of these gene targets were carried out. In brief, the raw data of 105 laryngeal cancer samples were firstly extracted and then standardized by log2 (count+1). After that, the optimal cutoff was determined using survminer package in R language based on the expressions of candidate genes. Subsequently, the samples were divided into high-and low-expression group according to the optimal cutoff. Following this, the survival analyses of gene targets were respectively preformed with survival package in R language. P < 0.05 was considered significant.

Predicting the binding mode between active compounds and key protein targets of BYQH decoction
The Protein Data Bank (PDB) offers the experimental data for determined three-dimensional (3D) biological macromolecules structure [26]. Herein, we would further explore the underlying molecular mechanism how active compounds interacted with target proteins of BYQH decoction. The PDB ID for the key protein was firstly acquired from PDB. Subsequently, the binding mode between key protein and corresponding active chemical ingredient was predicted by the LeDock tool. Finally, we used the Pymol software to visualize the predicted binding modes of key active ingredients to their corresponding targets.

Ethics approval and consent to participate
This study was approved by Ethics Committee of Provincial Hospital Affiliated to Shandong First Medical University and Shandong Provincial Hospital affiliated to Shandong University.

Availability of data and materials
The raw data supporting the conclusions of this manuscript will be made available by the AGING authors, without undue reservation, to any qualified researcher.

Highlights
The therapeutic mechanism of BYQH decoction on LSCC is analyzed by network pharmacology. Sitosterol and quercetin were key active ingredients BYQH decoction. TP53, EGFR, NOS3 and IL1B may be therapeutic targets of LSCC.

Target PPI network analysis
The potential interactions among target proteins of BYQH decoction was constructed and visualized by AGING STRING database. As shown in Figure 3, there were 128 protein nodes and 1282 edges in PPI network. The protein nodes with higher degree may have closer biological connection with other nodes. Herein, the top 15 nodes were regarded as hub genes and listed in Table 2

Cross-validation of target proteins
To further narrow the range of potential target genes of BYQH decoction, the candidate target genes were cross-validated. The genes related to LSCC were respectively acquired from Phenopedia and DisGeNET database, which provided information on gene-disease interactions. Then, the intersecting target genes were obtained between predicted genes and gene targets of BYQH decoction. Finally, a total of 10 intersectionassociated genes were extracted, including PTGS2, NOS3, BCL2, ADH1C, TP53, MPO, EGFR, GSTP1, IL1B and GSTM1 (Figure 4). Among them, five genes (TP53, EGFR, PTGS2, NOS3 and IL1B) served as the hub genes in PPI network. Therefore, these genes were considered to be key targets of BYQH decoction against LSCC.

Survival analysis
To explore the underlying impacts of 10 key gene targets on the clinical prognosis of LSCC, we conducted     Figure 5 and Supplementary Figure 1). Survival analysis revealed that seven genes (ADH1C, EGFR, GSTM1, GSTP1, IL1B, NOS3 and TP53) were significantly associated with the prognosis of HNSCC patients (P < 0.05; Figure 5). Moreover, HNSCC patients with a high expression level of ADH1C had a better prognosis (P = 0.012; Figure 5). However, the high expression levels of the remaining six genes were closely related to the unfavorable prognosis of HNSCC patients (EGFR: P = 0.014; GSTM1: P = 0.004; GSTP1: P = 0.034; IL1B: P = 0.03; NOS3: P = 0.039 and TP53: P = 0.02; Figure 5). Notably, TP53, EGFR, NOS3 and IL1B also acted as hub genes in PPI network.

Prediction of the binding mode between key bioactive compounds and proteins
The binding characteristics of key active ingredients of BYQH decoction and four key protein targets (TP53, EGFR, NOS3 and IL1B) were investigated. Herein, we focused on analyzing the binding modes between two overlapped active compounds (sitosterol and quercetin) among three Chinese herbs and four target proteins. Our analysis indicated that the corresponding protein targets of sitosterol were PGR, NCOA2, NR3C2, which were not targets studied. For quercetin, 75 target proteins were obtained, including four key protein targets. Subsequently, the chemical structure formula of quercetin was acquired from TSMCP database, and four target proteins were also searched from PDB database, including TP53 (PDB ID: 3ZME), EGFR (PDB ID: 1XKK), NOS3 (PDB ID: 6PP1) and IL1B (PDB ID: 5R86). Finally, the molecular docking was carried out using LeDock tool. Notably, a closer binding between proteins and small biological molecules indicated more energy released and a lower ΔG value. Figure 6 depicted the optimal binding modes of quercetin and its four target proteins.

DISCUSSION
TCM is an accepted medical practice and has been used for the treatment of complicated diseases such as cancers over past few decades in China [27,28]. BYQH decoction, a TCM prescription, has been clinically proven to be effective as adjuvant treatment against LSCC. In this study, we performed a network pharmacological analysis to explore pharmacodynamic effects and therapeutic mechanism of BYQH decoction for LSCC. We identified 41 key active ingredients of BYQH decoction. Of these, sitosterol and quercetin were intersected among three Chinese herbs of BYQH decoction. Moreover, these crucial components were corresponded to 137 target proteins. Top  Quercetin (3,3′,4′,5,7-pentahydroxyflavone), a common flavonol, is widely distributed in plant species, such as vegetables and grains. Our results showed that quercetin was a key active chemical compound of BYQH decoction. Extensive evidence suggests that quercetin exerts diverse biological functions, such as antioxidant, anti-inflammatory, anti-bacterial, antiviral and anticarcinogenic properties [29][30][31]. Numerous studies have evaluated the anticarcinogenic effects of quercetin and indicated that quercetin was involved in regulating several cancer-related pathways, such PI3K/Akt/ signaling pathways and MAPK/ERK1/2 pathways [32,33]. Sharma et al. report that quercetin induced human laryngeal HeP2 cells death and synergistically enhanced the antiproliferative ability of cisplatin in these cells [34]. Numerous researchers demonstrate that quercetin (50 μM) can significantly increase photodynamic therapy-induced cytotoxicity via reducing the cell viability of human larynx carcinoma cells (HEp-2) [35]. In addition, the apoptotic events are increased stronger in HEp-2 cells with the combination of quercetin than the single drug treatment [36]. These evidences reveals that quercetin may be a promising therapeutic molecule for LSCC treatment.
TP53, EGFR, NOS3 and IL1B were key hub genes in PPI network of target proteins of BYQH decoction and cross-validated in two databases (Phenopedia and DisGeNET). More notably, these four proteins were also target proteins of quercetin. p53 protein encoded by TP53 gene functions as a tumor-suppressive factor by regulating multiple cellular processes such as cell-cycle progression and apoptosis [37]. Overwhelming evidence suggests several mutations in TP53 can influence the progression of HNSCC and clinical treatment response [38,39]. Clemente-Soto AF et al. point out that quercetin can induce G2 phase cell cycle arrest and apoptosis in human cervical cancer cells, accompanied by upregulating p53 level [40]. Our molecular docking analysis suggested that quercetin can strongly bind with TP53. EGFR is frequently reported to be associated with HNSCC [41]. Chan et al. indicate that quercetin may suppress cell migration and invasion of HNSCC cells overexpressing EGFR by down-regulating the expression of MMP-2 and MMP-9 [42]. Similarly, Chung et al. indicate that there is a high frequency of EGFR copy number in HNSCC, which acts as a poor predictor for the prognosis of HNSCC patients [43]. IL1B is an inflammatory cytokine gene and acts as a therapeutic target for the treatment of head and neck cancer [44,45]. NOS3 is located on chromosome 7 (7q36) and regulates at transcriptional and posttranscriptional levels [46]. Previous studies identified numerous polymorphic sites of NOS3 such as single nucleotide polymorphism and insertion or deletion [47,48]. Moreover, the NOS3 polymorphisms play crucial roles in the molecular mechanism of cancers and clinical survivals of patients undergoing cancers, including laryngeal cancer [49][50][51]. Guo et al. recently AGING report that quercetin has a binding interaction with NOS3, which was consistent with our finding [52]. Besides, we also found the elevated level of NOS3 represented an unfavorable prognosis.
Besides, in this study, GO enrichment analysis showed that the four genes of interest (TP53, EGFR, NOS3 and IL1B) were significantly enriched in cellular response to chemical stimulus, response to chemical, response to AGING drug related biological processes. TP53 is a short-lived protein that plays a central role in mediating cellular response to stressful and genotoxic stimuli, such as anticancer drugs exposure [53]. EGFR has been proposed to be the prognostic marker for HNSCC that closely related with radiation sensitivity, tumor size and recurrence [54]. Additionally, our survival data showed that the high expression of TP53, EGFR, NOS3 and IL1B in HNSCC patients exhibited a poor prognosis, which was consistent with the previous reports mentioned above. Thus, we speculated that the quercetin might reduce drug resistance of HNSCC by targeting TP53, EGFR, NOS3 and IL1B.
Furthermore, PI3K-AKT signaling is considered to be a classical pathway related with various cancers [55,56]. The descending activation of PI3K-AKT signaling is closely associated with radiosensitivity enhancement of LSCC patients [57]. The reduced protein expressions of PI3K and AKT are accompanied by suppressed LSCC tumor growth [58]. In this study, PI3K-AKT signaling was identified to a significant pathway involved with TP53, EGFR, NOS3. Our findings supported the significant role of PI3K-AKT signaling pathway in LSCC and proposed the target role of the key proteins identified in this study.
Although our study identified several key components and target proteins of BYQH decoction, relevant experimental assays are needed to validate our findings. Although were considered as the prognostic indicators for patients with LSCC, The expression patterns of the four candidate prognostic indicators in LSCC have not been investigated in clinical samples. Moreover, the relationships of these prognostic genes and some clinical parameters are also not analyzed. Additionally, many chemical compounds are possibly not identified because the limited data is available in public databases, which may lead to some false positive results.

CONCLUSIONS
Network pharmacology analysis was carried out to elucidate the therapeutic role of BYQH decoction on LSCC. Quercetin was predicted to be the active compounds of BYQH decoction by targeting TP53, EGFR, NOS3 and IL1B involved in drug resistance and PI3K-AKT signaling pathway. Moreover, the four genes of TP53, EGFR, NOS3 and IL1B exerted promising prognostic potential for LSCC. Our results offered an important clue for uncovering the underlying mechanism of BYQH decoction in the treatment of LSCC. However, more experimental evidences about the effects of quercetin on LSCC in vivo and in vitro by targeting TP53, EGFR, NOS3 and IL1B are warranted in the near future.

AUTHOR CONTRIBUTIONS
Aiai Lv, Xuanchen Zhou and Kun Gao carried out the Conception and design of the research, Kun Gao and Xianwei Gong participated in the Acquisition of data. Kun Gao and Yanan Zhu carried out the Analysis and interpretation of data. Kun Gao and Hui Wang participated in the design of the study and performed the statistical analysis. Zhiyong Yue participated in Obtaining funding. Aiai Lv, Xuanchen Zhou, Kun Gao and Aiai Lv conceived of the study, and participated in its design and coordination and helped to draft the manuscript and revision of manuscript for important intellectual content. All authors read and approved the final manuscript.

CONFLICTS OF INTEREST
The authors declare no conflicts of interest related to this study.