# Investigating the Mechanism ofScutellariae barbata Herba in the Treatment of Colorectal Cancer by Network Pharmacology and Molecular Docking

Authors: Xiangjun Qi; Hongbin Xu; Peng Zhang; Guoming Chen; Zhiqiang Chen; Caishan Fang; Lizhu Lin
Journal: Evidence-Based Complementary and Alternative Medicine
 (2021)

---

## Abstract

Background. Colorectal cancer (CRC) is one of the most common gastrointestinal tumors, which accounts for approximately 10% of all diagnosed cancers and cancer deaths worldwide per year. Scutellariae barbatae Herba (SBH) is one of the most frequently used traditional Chinese medicine (TCM) in the treatment of CRC. Although many experiments have been carried out to explain the mechanisms of SBH, the mechanisms of SBH have not been illuminated fully. Thus, we constructed a network pharmacology and molecular docking to investigate the mechanisms of SBH. Methods. We adopted active constituent prescreening, target predicting, protein-protein interaction (PPI) analysis, Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, differentially expressed gene analysis, and molecular docking to establish a system pharmacology database of SBH against CRC. Results. A total of 64 active constituents of SBH were obtained and 377 targets were predicted, and the result indicated that quercetin, luteolin, wogonin, and apigenin were the main active constituents of SBH. Glucocorticoid receptor (NR3C1), pPhosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit alpha isoform (PIK3CA), cellular tumor antigen p53 (TP53), transcription factor AP-1 (JUN), mitogen-activated protein kinase 1 (MAPK1), Myc protooncogene protein (MYC), cyclin-dependent kinase 1 (CDK1), and broad substrate specificity ATP-binding cassette transporter ABCG2 (ABCG2) were the major targets of SBH in the treatment of CRC. GO analysis illustrated that the core biological process regulated by SBH was the regulation of the cell cycle. Thirty pathways were presented and 8 pathways related to CRC were involved. Molecular docking presented the binding details of 3 key targets with 6 active constituents. Conclusions. The mechanisms of SBH against CRC depend on the synergistic effect of multiple active constituents, multiple targets, and multiple pathways.


---

## Body

## 1. Introduction

Colorectal cancer (CRC) is one of the most common gastrointestinal tumors. Nowadays, CRC accounts for approximately 10% of all diagnosed cancers and cancer deaths worldwide per year [1]. The amounts of CRC patients will increase to 2.5 million in 2035 around the world according to the prediction of WHO [1, 2]. Previous studies have indicated that the majority of CRC cells originate from colon stem cells, which were located at the base of colonic crypts [3, 4], and several risk factors such as smoking [5], excessive alcohol ingestion [6], and red meat consumption [7] have been identified to contribute to the genetic mutations of colon stem cells. The therapies for CRC include endoscopic treatment, surgery, radiotherapy, chemotherapy, and immunotherapy [1]. Although there is brilliant progress in diagnostic techniques and treatment strategies [8], there still exist a sizeable percentage of patients at an advanced stage of CRC along with a high degree of metastasis when diagnosed [9]. In addition, the cost of anticancer drugs keeps increasing, which brings heavy economic burden to families, especially in developing countries [10]. Thus, it is necessary to search for more cost-effective and less toxic drugs.Traditional Chinese medicine (TCM) has a 5000-year history in China, and it is presently regarded as a useful complementary and alternative medication worldwide [11]. Many historic studies have shown that TCM could induce apoptosis, tumor growth, angiogenesis, and metastasis [12–15]. TCM also plays a role in ameliorating the side effects engendered by radiotherapy and chemotherapy [16]. Scutellariae barbatae Herba (also known as BanZhiLian in Chinese) is the dried full plant of Scutellaria barbata D. Don in the Lamiaceae family, which has been used for thousands of years in China as a “heat-clearing and detoxifying” drug [17, 18]. Moreover, the anticancer effect of SBH has been examined in many types of cancer including lung cancer, hepatic cancer, breast cancer, colorectal cancer, leukemia, and prostate cancer [17]. Although a considerable amount of literature has been published on searching the mechanisms of SBH, the potential mechanisms of SBH have not yet been systematically investigated.The traditional drug discovery paradigm focuses on screening exquisite ligands, and the presupposition of “one drug for one target for one disease” is accepted by the majority of researchers [19]. However, this paradigm is not suitable for TCM, which is characterized by multiple components, multiple targets, and multiple pathways [20]. Network pharmacology highlights the consolidation of drug targets, biological network, and pharmacology network [19], which provides a feasible method for TCM exploration. It has already been used for searching core active constituents and targets and detecting the mechanisms of TCM.Molecular docking is an essential tool in computer-assisted drug design, which was firstly proposed in the mid-80s with the purpose of predicting mode of ligand and protein and virtually screening digital compound libraries to reduce expense and speed up drug discovery. The accuracy of molecular docking keeps increasing with the development of computing power and hardware capability.In this study, we adopted network pharmacology and molecular docking to investigate the mechanisms of SBH in the treatment of CRC. Our work involved 6 parts as follows: (1) collecting active constituents of SBH from online TCM databases and literatures; (2) predicting the targets of SBH and disease targets of CRC; (3) constructing network to show the interaction between SBH and CRC; (4) conducting function analysis including Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis; (5) seeking differentially expressed genes of CRC with a bioinformatics analysis; (6) implementing molecular docking after a literature retrieval.


## 2. Materials and Methods

A detailed diagram used to describe the overall design of this study is shown in Figure1.Figure 1
Overall design of the study.


### 2.1. Data Sources of Network Pharmacology

#### 2.1.1. Active Constituents of SBH and Target Proteins Prediction

The active constituents of SBH were obtained in 3 steps. Step 1: data retrieval was performed with 7 databases: Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP,http://lsp.nwu.edu.cn/tcmsp.php) [21], The Encyclopedia of Traditional Chinese Medicine (ETCM, http://www.nrc.ac.cn:9090/ETCM/) [22], Traditional Chinese Medicines Integrated Database (TCMID, http://www.megabionet.org/tcmid/) [23], Traditional Chinese Medicine Information Database (TCM-ID, http://bidd.nus.edu.sg/group/TCMsite/) [24], TCM-Mesh System (http://mesh.tcm.microbioinformatics.org/) [25], TCMGeneDIT Databases (http://tcm.lifescience.ntu.edu.tw/index.html) [26], and Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine (BATMAN-TCM, http://bionet.ncpsb.org/batman-tcm/) [27]. The benefit of this approach is that all the databases provide a comprehensive study for Chinese herb. The oral bioavailability (OB) and drug-likeness (DL) were adopted to find out the active constituents with better pharmacokinetics [28], and active constituents with OB ≥ 30% and DL ≥ 0.18 were selected for subsequent network construction [29]. Step 2: Anticancer Herbs Database of Systems Pharmacology (CancerHSP, http://lsp.nwsuaf.edu.cn/CancerHSP.php) [30] was retrieved to collect the active constituents with anticancer activity, which contains 2439 anticancer herbal medicines with 3575 anticancer active constituents. In addition, it also provides 832 targets of the active constituents that are predicted by state-of-art methods or collected from the literature. Step 3: to obtain the active constituents which showed bioavailability activities for CRC but were not recorded in aforementioned databases, we constructed a text mining in PubMed using “Scutellariae Barbata Herba” and “cancer” as search terms. After manual filtering of the search results, barbatin F, barbatin G [31], barbatin H [32], SPS2p [33], and SBPW3 [34] were supplemented.Target protein prediction was based on TCMSP database and DrugBank database (https://www.drugbank.ca/) [35]. Active constituents, whose targets could not be predicted in TCMSP or DrugBank, were predicted by SwissTargetPrediction (http://www.swisstargetprediction.ch), and the top 15 predicted targets were selected for each result of the predicted target classes [36, 37].To better visualize the relationship between active constituents and target proteins, we coded each active constituent and the SBH active constituent-target protein network was generated using Cytoscape (v3.7.1) software (http://www.cytoscape.org/) [38]. We used Network Analyzer (a cytoscape plugin) to calculate the degree and betweenness centrality of the network. The degree of an active constituent n is the number of target proteins linked to n. The betweenness centrality reflects the extent to which a node acts as a communication intermediate in the network.


#### 2.1.2. CRC-Related Genes and Corresponding Proteins

CRC-related genes were obtained from 4 databases: GeneCards: the Human Gene Database (https://www.genecards.org/), which integrates more than 190 data source about gene, disease, pathway, and compound [39]; Online Mendelian Inheritance in Man (OMIM, https://omim.org/), which contains information on all known Mendelian disorders and over 15000 genes [40]; Therapeutic Target Database (TTD, http://db.idrblab.net/ttd/), which is a database to provide therapeutic proteins and nucleic acid targets, the targeted disease, and the correspondence between drugs and these targets [41]; Comparative Toxicogenomics Database (CTD, http://ctdbase.org/), which supplies information of 3 interactions: chemical-gene/protein, chemical-disease, and gene-disease [42]. Four databases were searched with the same keywords: “colorectal cancer” or “colorectal Neoplasms” or “Colorectal Carcinoma” or “Colorectal Tumor”, and the species was limited to Homo sapiens. The proteins corresponding to the CRC-related genes were standardized using UniProt database for subsequent analysis.


#### 2.1.3. Potential Active Target Proteins (PATPs)

Finally, we took the intersection between the target proteins of SBH and CRC-related proteins as the PATPs for subsequent analysis.


### 2.2. Data Analysis of Network Pharmacology

#### 2.2.1. Protein-Protein Interaction (PPI) Analysis

To figure out the interactions between PATPs, we used the STRING database (https://string-db.org/) to construct the PPI network [43], the species was limited to Homo sapiens, and we selected confidence score as 0.95 in the minimum (low confidence < 0.4, medium ≤ 0.7, and high > 0.7). The PPI data was exported as a.tsv file for further analysis. Afterward, we used CytoNCA (a Cytoscape plugin) to evaluate the PPI network [44], and the PPI network was constructed by the top 150 proteins.


#### 2.2.2. Gene Ontology (GO) Analysis

The GO enrichment analysis is one of the most common procedures for determining potential molecular mechanisms of drugs. We did GO biological process (BP) analysis using MCODE (a Cytoscape plugin) [45] and ClueGO (a Cytoscape plugin) [46]. First, the PPI data were filtered in Cytoscape using MCODE plugin and the top 3 clusters (sorted by cluster score) were analyzed by the ClueGo plugin. During this procedure, the significance level of GO terms was set to 0.05, and the species was limited to Homo sapiens.


#### 2.2.3. Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis

Traditionally, one of the most well-known methods of exploring biological pathways and potential biological functions is KEGG pathway enrichment analysis. In our study, KEGG pathway analysis was performed withClusterProfiler package of R language [47], and the significance level of GO terms was set to 0.05 and species was limited to Homo sapiens.


### 2.3. Bioinformatics Analysis

#### 2.3.1. Microarray Data

Microarray data was collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) with the following limits: “colorectal cancer” (keyword), “Homo sapiens” (study organism), and “tissue” (attribute name). After collection, the GSE32323 dataset, which contained 17 pairs of matched CRC samples based on the GPL570 platform, and the GSE44076 dataset, which contained 98 pairs of matched CRC samples based on the GPL13667 platform, were selected for seeking differentially expressed genes of CRC.


#### 2.3.2. Identification of Differentially Expressed Genes

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), an online web tool based on limma package of R language, is one of the most common tools of determining differentially expressed genes (DEGs) of upregulation or downregulation. In our study, we used GEO2R to finish this procedure, and the results were presented as a table of genes. We set adjusted p value <0.05 and log2fold changeFC≥2 as a criterion to screen out DEGs. The proteins corresponding to the DEGs are standardized using UniProt database for subsequent molecular docking.


### 2.4. Molecular Docking between Key Targets and Active Constituents

#### 2.4.1. Determination of Key Target Proteins and Active Constituents

We took the intersection between PATPs and DEGs-related proteins as the key target proteins for molecular docking. Afterward, we plotted a network diagram to represent the correspondence between key target proteins and active constituents of SBH. Finally, a literature review was conducted to learn the relationship between key target proteins and CRC, and those that were experimentally validated as potential therapeutic targets for CRC would be selected for molecular docking.


#### 2.4.2. Molecular Docking Simulation

First, a suitable protein structure of each target was obtained from the RCSB Protein Data Bank (https://www.rcsb.org/) [48], and the suitable protein was required to satisfy the following 3 conditions as far as possible: (1) the protein owned a 3D structure with a high resolution; (2) the protein owned one or more original ligands; (3) the ligand owned a similar structure with active constituent. Second, Chimera (v1.14) software was applied to remove the heteroatom and water molecule from the proteins and divide the proteins into ligands and receptors. The ligand and receptor files were converted into a pdbqt file then by AutoDockTools (v1.5.6) software. Third, the 2D structures of active constituents were downloaded from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/) [49]. AutoDockTools (v1.5.6) software was used to convert them into a pdbqt file then. Fourth, a grid box size set as 40 × 40 × 40 points with a Vina spacing of 1.0 Å was generated. Finally, we used AutoDock Vina (v1.1.2), an open-source program for molecular docking simulation, to accomplish the docking stimulation, which significantly improved the average accuracy of the binding mode predictions and speed of docking when compared with AutoDock [50]. The result of docking was analyzed by Discovery Studio.


## 2.1. Data Sources of Network Pharmacology

### 2.1.1. Active Constituents of SBH and Target Proteins Prediction

The active constituents of SBH were obtained in 3 steps. Step 1: data retrieval was performed with 7 databases: Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP,http://lsp.nwu.edu.cn/tcmsp.php) [21], The Encyclopedia of Traditional Chinese Medicine (ETCM, http://www.nrc.ac.cn:9090/ETCM/) [22], Traditional Chinese Medicines Integrated Database (TCMID, http://www.megabionet.org/tcmid/) [23], Traditional Chinese Medicine Information Database (TCM-ID, http://bidd.nus.edu.sg/group/TCMsite/) [24], TCM-Mesh System (http://mesh.tcm.microbioinformatics.org/) [25], TCMGeneDIT Databases (http://tcm.lifescience.ntu.edu.tw/index.html) [26], and Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine (BATMAN-TCM, http://bionet.ncpsb.org/batman-tcm/) [27]. The benefit of this approach is that all the databases provide a comprehensive study for Chinese herb. The oral bioavailability (OB) and drug-likeness (DL) were adopted to find out the active constituents with better pharmacokinetics [28], and active constituents with OB ≥ 30% and DL ≥ 0.18 were selected for subsequent network construction [29]. Step 2: Anticancer Herbs Database of Systems Pharmacology (CancerHSP, http://lsp.nwsuaf.edu.cn/CancerHSP.php) [30] was retrieved to collect the active constituents with anticancer activity, which contains 2439 anticancer herbal medicines with 3575 anticancer active constituents. In addition, it also provides 832 targets of the active constituents that are predicted by state-of-art methods or collected from the literature. Step 3: to obtain the active constituents which showed bioavailability activities for CRC but were not recorded in aforementioned databases, we constructed a text mining in PubMed using “Scutellariae Barbata Herba” and “cancer” as search terms. After manual filtering of the search results, barbatin F, barbatin G [31], barbatin H [32], SPS2p [33], and SBPW3 [34] were supplemented.Target protein prediction was based on TCMSP database and DrugBank database (https://www.drugbank.ca/) [35]. Active constituents, whose targets could not be predicted in TCMSP or DrugBank, were predicted by SwissTargetPrediction (http://www.swisstargetprediction.ch), and the top 15 predicted targets were selected for each result of the predicted target classes [36, 37].To better visualize the relationship between active constituents and target proteins, we coded each active constituent and the SBH active constituent-target protein network was generated using Cytoscape (v3.7.1) software (http://www.cytoscape.org/) [38]. We used Network Analyzer (a cytoscape plugin) to calculate the degree and betweenness centrality of the network. The degree of an active constituent n is the number of target proteins linked to n. The betweenness centrality reflects the extent to which a node acts as a communication intermediate in the network.


### 2.1.2. CRC-Related Genes and Corresponding Proteins

CRC-related genes were obtained from 4 databases: GeneCards: the Human Gene Database (https://www.genecards.org/), which integrates more than 190 data source about gene, disease, pathway, and compound [39]; Online Mendelian Inheritance in Man (OMIM, https://omim.org/), which contains information on all known Mendelian disorders and over 15000 genes [40]; Therapeutic Target Database (TTD, http://db.idrblab.net/ttd/), which is a database to provide therapeutic proteins and nucleic acid targets, the targeted disease, and the correspondence between drugs and these targets [41]; Comparative Toxicogenomics Database (CTD, http://ctdbase.org/), which supplies information of 3 interactions: chemical-gene/protein, chemical-disease, and gene-disease [42]. Four databases were searched with the same keywords: “colorectal cancer” or “colorectal Neoplasms” or “Colorectal Carcinoma” or “Colorectal Tumor”, and the species was limited to Homo sapiens. The proteins corresponding to the CRC-related genes were standardized using UniProt database for subsequent analysis.


### 2.1.3. Potential Active Target Proteins (PATPs)

Finally, we took the intersection between the target proteins of SBH and CRC-related proteins as the PATPs for subsequent analysis.


## 2.1.1. Active Constituents of SBH and Target Proteins Prediction

The active constituents of SBH were obtained in 3 steps. Step 1: data retrieval was performed with 7 databases: Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP,http://lsp.nwu.edu.cn/tcmsp.php) [21], The Encyclopedia of Traditional Chinese Medicine (ETCM, http://www.nrc.ac.cn:9090/ETCM/) [22], Traditional Chinese Medicines Integrated Database (TCMID, http://www.megabionet.org/tcmid/) [23], Traditional Chinese Medicine Information Database (TCM-ID, http://bidd.nus.edu.sg/group/TCMsite/) [24], TCM-Mesh System (http://mesh.tcm.microbioinformatics.org/) [25], TCMGeneDIT Databases (http://tcm.lifescience.ntu.edu.tw/index.html) [26], and Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine (BATMAN-TCM, http://bionet.ncpsb.org/batman-tcm/) [27]. The benefit of this approach is that all the databases provide a comprehensive study for Chinese herb. The oral bioavailability (OB) and drug-likeness (DL) were adopted to find out the active constituents with better pharmacokinetics [28], and active constituents with OB ≥ 30% and DL ≥ 0.18 were selected for subsequent network construction [29]. Step 2: Anticancer Herbs Database of Systems Pharmacology (CancerHSP, http://lsp.nwsuaf.edu.cn/CancerHSP.php) [30] was retrieved to collect the active constituents with anticancer activity, which contains 2439 anticancer herbal medicines with 3575 anticancer active constituents. In addition, it also provides 832 targets of the active constituents that are predicted by state-of-art methods or collected from the literature. Step 3: to obtain the active constituents which showed bioavailability activities for CRC but were not recorded in aforementioned databases, we constructed a text mining in PubMed using “Scutellariae Barbata Herba” and “cancer” as search terms. After manual filtering of the search results, barbatin F, barbatin G [31], barbatin H [32], SPS2p [33], and SBPW3 [34] were supplemented.Target protein prediction was based on TCMSP database and DrugBank database (https://www.drugbank.ca/) [35]. Active constituents, whose targets could not be predicted in TCMSP or DrugBank, were predicted by SwissTargetPrediction (http://www.swisstargetprediction.ch), and the top 15 predicted targets were selected for each result of the predicted target classes [36, 37].To better visualize the relationship between active constituents and target proteins, we coded each active constituent and the SBH active constituent-target protein network was generated using Cytoscape (v3.7.1) software (http://www.cytoscape.org/) [38]. We used Network Analyzer (a cytoscape plugin) to calculate the degree and betweenness centrality of the network. The degree of an active constituent n is the number of target proteins linked to n. The betweenness centrality reflects the extent to which a node acts as a communication intermediate in the network.


## 2.1.2. CRC-Related Genes and Corresponding Proteins

CRC-related genes were obtained from 4 databases: GeneCards: the Human Gene Database (https://www.genecards.org/), which integrates more than 190 data source about gene, disease, pathway, and compound [39]; Online Mendelian Inheritance in Man (OMIM, https://omim.org/), which contains information on all known Mendelian disorders and over 15000 genes [40]; Therapeutic Target Database (TTD, http://db.idrblab.net/ttd/), which is a database to provide therapeutic proteins and nucleic acid targets, the targeted disease, and the correspondence between drugs and these targets [41]; Comparative Toxicogenomics Database (CTD, http://ctdbase.org/), which supplies information of 3 interactions: chemical-gene/protein, chemical-disease, and gene-disease [42]. Four databases were searched with the same keywords: “colorectal cancer” or “colorectal Neoplasms” or “Colorectal Carcinoma” or “Colorectal Tumor”, and the species was limited to Homo sapiens. The proteins corresponding to the CRC-related genes were standardized using UniProt database for subsequent analysis.


## 2.1.3. Potential Active Target Proteins (PATPs)

Finally, we took the intersection between the target proteins of SBH and CRC-related proteins as the PATPs for subsequent analysis.


## 2.2. Data Analysis of Network Pharmacology

### 2.2.1. Protein-Protein Interaction (PPI) Analysis

To figure out the interactions between PATPs, we used the STRING database (https://string-db.org/) to construct the PPI network [43], the species was limited to Homo sapiens, and we selected confidence score as 0.95 in the minimum (low confidence < 0.4, medium ≤ 0.7, and high > 0.7). The PPI data was exported as a.tsv file for further analysis. Afterward, we used CytoNCA (a Cytoscape plugin) to evaluate the PPI network [44], and the PPI network was constructed by the top 150 proteins.


### 2.2.2. Gene Ontology (GO) Analysis

The GO enrichment analysis is one of the most common procedures for determining potential molecular mechanisms of drugs. We did GO biological process (BP) analysis using MCODE (a Cytoscape plugin) [45] and ClueGO (a Cytoscape plugin) [46]. First, the PPI data were filtered in Cytoscape using MCODE plugin and the top 3 clusters (sorted by cluster score) were analyzed by the ClueGo plugin. During this procedure, the significance level of GO terms was set to 0.05, and the species was limited to Homo sapiens.


### 2.2.3. Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis

Traditionally, one of the most well-known methods of exploring biological pathways and potential biological functions is KEGG pathway enrichment analysis. In our study, KEGG pathway analysis was performed withClusterProfiler package of R language [47], and the significance level of GO terms was set to 0.05 and species was limited to Homo sapiens.


## 2.2.1. Protein-Protein Interaction (PPI) Analysis

To figure out the interactions between PATPs, we used the STRING database (https://string-db.org/) to construct the PPI network [43], the species was limited to Homo sapiens, and we selected confidence score as 0.95 in the minimum (low confidence < 0.4, medium ≤ 0.7, and high > 0.7). The PPI data was exported as a.tsv file for further analysis. Afterward, we used CytoNCA (a Cytoscape plugin) to evaluate the PPI network [44], and the PPI network was constructed by the top 150 proteins.


## 2.2.2. Gene Ontology (GO) Analysis

The GO enrichment analysis is one of the most common procedures for determining potential molecular mechanisms of drugs. We did GO biological process (BP) analysis using MCODE (a Cytoscape plugin) [45] and ClueGO (a Cytoscape plugin) [46]. First, the PPI data were filtered in Cytoscape using MCODE plugin and the top 3 clusters (sorted by cluster score) were analyzed by the ClueGo plugin. During this procedure, the significance level of GO terms was set to 0.05, and the species was limited to Homo sapiens.


## 2.2.3. Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis

Traditionally, one of the most well-known methods of exploring biological pathways and potential biological functions is KEGG pathway enrichment analysis. In our study, KEGG pathway analysis was performed withClusterProfiler package of R language [47], and the significance level of GO terms was set to 0.05 and species was limited to Homo sapiens.


## 2.3. Bioinformatics Analysis

### 2.3.1. Microarray Data

Microarray data was collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) with the following limits: “colorectal cancer” (keyword), “Homo sapiens” (study organism), and “tissue” (attribute name). After collection, the GSE32323 dataset, which contained 17 pairs of matched CRC samples based on the GPL570 platform, and the GSE44076 dataset, which contained 98 pairs of matched CRC samples based on the GPL13667 platform, were selected for seeking differentially expressed genes of CRC.


### 2.3.2. Identification of Differentially Expressed Genes

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), an online web tool based on limma package of R language, is one of the most common tools of determining differentially expressed genes (DEGs) of upregulation or downregulation. In our study, we used GEO2R to finish this procedure, and the results were presented as a table of genes. We set adjusted p value <0.05 and log2fold changeFC≥2 as a criterion to screen out DEGs. The proteins corresponding to the DEGs are standardized using UniProt database for subsequent molecular docking.


## 2.3.1. Microarray Data

Microarray data was collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) with the following limits: “colorectal cancer” (keyword), “Homo sapiens” (study organism), and “tissue” (attribute name). After collection, the GSE32323 dataset, which contained 17 pairs of matched CRC samples based on the GPL570 platform, and the GSE44076 dataset, which contained 98 pairs of matched CRC samples based on the GPL13667 platform, were selected for seeking differentially expressed genes of CRC.


## 2.3.2. Identification of Differentially Expressed Genes

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/), an online web tool based on limma package of R language, is one of the most common tools of determining differentially expressed genes (DEGs) of upregulation or downregulation. In our study, we used GEO2R to finish this procedure, and the results were presented as a table of genes. We set adjusted p value <0.05 and log2fold changeFC≥2 as a criterion to screen out DEGs. The proteins corresponding to the DEGs are standardized using UniProt database for subsequent molecular docking.


## 2.4. Molecular Docking between Key Targets and Active Constituents

### 2.4.1. Determination of Key Target Proteins and Active Constituents

We took the intersection between PATPs and DEGs-related proteins as the key target proteins for molecular docking. Afterward, we plotted a network diagram to represent the correspondence between key target proteins and active constituents of SBH. Finally, a literature review was conducted to learn the relationship between key target proteins and CRC, and those that were experimentally validated as potential therapeutic targets for CRC would be selected for molecular docking.


### 2.4.2. Molecular Docking Simulation

First, a suitable protein structure of each target was obtained from the RCSB Protein Data Bank (https://www.rcsb.org/) [48], and the suitable protein was required to satisfy the following 3 conditions as far as possible: (1) the protein owned a 3D structure with a high resolution; (2) the protein owned one or more original ligands; (3) the ligand owned a similar structure with active constituent. Second, Chimera (v1.14) software was applied to remove the heteroatom and water molecule from the proteins and divide the proteins into ligands and receptors. The ligand and receptor files were converted into a pdbqt file then by AutoDockTools (v1.5.6) software. Third, the 2D structures of active constituents were downloaded from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/) [49]. AutoDockTools (v1.5.6) software was used to convert them into a pdbqt file then. Fourth, a grid box size set as 40 × 40 × 40 points with a Vina spacing of 1.0 Å was generated. Finally, we used AutoDock Vina (v1.1.2), an open-source program for molecular docking simulation, to accomplish the docking stimulation, which significantly improved the average accuracy of the binding mode predictions and speed of docking when compared with AutoDock [50]. The result of docking was analyzed by Discovery Studio.


## 2.4.1. Determination of Key Target Proteins and Active Constituents

We took the intersection between PATPs and DEGs-related proteins as the key target proteins for molecular docking. Afterward, we plotted a network diagram to represent the correspondence between key target proteins and active constituents of SBH. Finally, a literature review was conducted to learn the relationship between key target proteins and CRC, and those that were experimentally validated as potential therapeutic targets for CRC would be selected for molecular docking.


## 2.4.2. Molecular Docking Simulation

First, a suitable protein structure of each target was obtained from the RCSB Protein Data Bank (https://www.rcsb.org/) [48], and the suitable protein was required to satisfy the following 3 conditions as far as possible: (1) the protein owned a 3D structure with a high resolution; (2) the protein owned one or more original ligands; (3) the ligand owned a similar structure with active constituent. Second, Chimera (v1.14) software was applied to remove the heteroatom and water molecule from the proteins and divide the proteins into ligands and receptors. The ligand and receptor files were converted into a pdbqt file then by AutoDockTools (v1.5.6) software. Third, the 2D structures of active constituents were downloaded from the PubChem website (https://pubchem.ncbi.nlm.nih.gov/) [49]. AutoDockTools (v1.5.6) software was used to convert them into a pdbqt file then. Fourth, a grid box size set as 40 × 40 × 40 points with a Vina spacing of 1.0 Å was generated. Finally, we used AutoDock Vina (v1.1.2), an open-source program for molecular docking simulation, to accomplish the docking stimulation, which significantly improved the average accuracy of the binding mode predictions and speed of docking when compared with AutoDock [50]. The result of docking was analyzed by Discovery Studio.


## 3. Results

### 3.1. SBH Active Constituent-Target Network

In our study, we obtained 64 active constituents in total. All active constituents were ranked and detailed information of each active constituent is shown in Table1. The result of SBH active constituent-target network is presented in Supplementary Figure S1, which was composed of 441 nodes (including 64 active constituent nodes and 377 target nodes) and 1161 edges, and the result of the top 10 active constituents and targets according to the degree is set out in Supplementary Table S1.Table 1
Active constituents ofScutellariae barbata herba.
RankingActive constituentsMol01(2R)-5,7-Dihydroxy-2-(4-hydroxyphenyl)chroman-4-oneMol0224-Ethylcholest-4-en-3-oneMol035,7,4′-Trihydroxy-6-methoxyflavanoneMol045,7,4′-Trihydroxy-8-methoxyflavanoneMol055-Hydroxy-7,8-dimethoxy-2-(4-methoxyphenyl)chromoneMol066-HydroxynaringeninMol077-Hydroxy-5,8-dimethoxy-2-phenyl-chromoneMol089,19-Cyclolanost-24-en-3-olMol09BaicaleinMol10BaicalinMol11Beta-sitosterolMol12CampesterolMol13CarthamidinMol14Chrysin-5-methyletherMol15CLRMol16Daucostero_qtMol17DinatinMol18EriodictyolMol19LuteolinMol20MoslosooflavoneMol21QuercetinMol22RhamnazinMol23RivularinMol24SalvigeninMol25SitosterolMol26Sitosterol acetateMol27Stigmasta-5,22-dien-3-ol-acetateMol28Stigmastan-3,5,22-trieneMol29StigmasterolMol30WogoninMol31Scutebarbatine FMol326-O-Nicotinoylscutebarbatine GMol33ApigeninMol34Scutelinquanine DMol35Scutebarbatine EMol36Scutebarbatine NMol376-O-(2-Carbonyl-3-methylbutanoyl)scutehenanine AMol386-O-Acetylscutehenanine AMol397-O-Nicotinoylscutebarbatine HMol40Barbatellarine BMol41Barbatin CMol42Barbatin EMol43Scutebarbatine CMol44Scutebarbatine DMol45Scutebarbatine HMol46Scutebata AMol47Scutehenanine AMol48Scutehenanine DMol496,7-Di-O-nicotinoylscutebarbatine GMol506-O-Nicotinoyl-7-O-acetylscutebarbatine GMol51Barbatin AMol52Barbatin BMol53Barbatin DMol54Pheophorbide aMol55Scutebarbatine BMol56Scutehenanine BMol57Scutehenanine CMol58Scutebarbatine OMol59Scutebarbatine GMol60ScutellarinMol617-Acetoxybarbatin CMol626-Acetoxybarbatin CMol63Barbatin FMol64Barbatin H


### 3.2. Potential Active Target Proteins and PPI Network

A total of 377 targets of SBH and 6390 target genes of CRC were obtained, and the intersection of the two groups was regarded as PATPs (Figure2). A total of 297 proteins were identified and used for constructing the PPI network. After limiting the species to Homo sapiens and setting a minimum confidence score as 0.95, 150 proteins were retained as shown in Figure 3. TP53, JUN, AKT1, MAPK1, PIK3CA, etc. were the core proteins of PPI.Figure 2
Potential active target proteins.Figure 3
PPI network of potential active target proteins. The nodes represent proteins, and the color of the node represents the degree value.


### 3.3. GO Enrichment Analysis

We constructed a cluster analysis using MCODE before GO enrichment analysis, and 9 clusters were obtained altogether (Table2). Afterwards, the top 3 clusters with high cluster scores were used for the procedure of GO enrichment analysis as shown in Figure 4. Interestingly, the majority of biological processes in cluster 1 were described as follows: negative regulation of G1/S transition of mitotic cell cycle (45.16%), cyclin-dependent protein serine/threonine kinase regulator activity (25.81%), and regulation of G1/S transition of mitotic cell cycle (19.35%), and they are of great importance to cell proliferation and influence the occurrence of cancer. Cluster 2 was constructed by membrane protein ectodomain proteolysis and notch receptor processing, which is ligand-dependent. Cluster 3 was composed of chemokine-mediated signaling pathway (40.51%), lipopolysaccharide-mediated signaling pathway (16.46%), mononuclear cell migration (12.66%), etc.Table 2
Detailed information of 9 clusters using MCODE.
ClusterScoreNodesEdges17.846145125.661434.211204043.4111753.3334563.33345733382.933162292.8891013Figure 4
GO enrichment analysis of clusters 1–3 (a–c). For the whole figure, different colors represent different biological process groups. The network diagrams represent the connections between biological process groups, the pie charts represent the proportion of each group, and the bar charts represent the detailed term of each group.
(a)(b)(c)


### 3.4. KEGG Pathway Analysis

KEGG pathway analysis was conducted with theClusterProfiler package of R language, and 174 significant pathways (adjusted p value < 0.05) were identified. We sorted them with adjusted p value and the top 30 terms are shown in Figure 5(a). Eight KEGG pathways related to human solid cancers including prostate cancer, small cell lung cancer, pancreatic cancer, bladder cancer, non-small cell lung cancer, hepatocellular carcinoma, colorectal cancer, and breast cancer were gathered. KEGG website summarized 9 pathways associated with CRC, eight of which were significantly enriched in our results (Figure 5(b)). Notably, the result of the colorectal cancer pathway is significant with adjusted p value < 0.05.Figure 5
KEGG pathway analysis of potential active targets. (a) For the top 30 terms sorted by adjustedp value and (b) for the eight pathways related to CRC summarized by KEGG website. The color of each term represents the different adjusted p values <  0.05 and the size of each term represents the gene counts.
(a)(b)


### 3.5. Identification of Key Target Protein and Active Constituent

In this section, we did an analysis of DEGs in CRC firstly using GEO2R. Two datasets including GSE32323 and GSE44076 were screened from the GEO database. A total of 369 DEGs were obtained from the GSE32323 dataset and 541 DEGs were obtained from the GSE44076. We normalized the proteins corresponding to these DEGs by the STRING database. Seventeen key target proteins were remaining after intersecting with PATPs (Figure6): ABCG2, SFRP1, MMP7, MYC, CA2, NR3C2, HSD17B2, PLAU, ADH1C, HSD11B2, MMP3, XDH, FABP6, CDK1, MMP1, MMP12, and SPP1. We constructed a network to show the active constituents that corresponded to these target proteins as presented in Figure 7. Among the active constituents mapped by key target proteins, quercetin, luteolin, baicalein, etc., were identified to be the top 10 active constituents in terms of degree value in the SBH active constituent-target network. After a literature review of these key target proteins, we chose MYC, ABCG2, and CDK1 for molecular docking.Figure 6
Intersections of differentially expressed genes related proteins and potential active target proteins.Figure 7
The network of key target proteins and related active constituents. The circular nodes represent the key target proteins, and the diamond nodes represent the active constituents. The color of the node represents the degree value.


### 3.6. Molecular Docking Analysis

To verify how an active constituent binds to target as previously referred to a molecular docking using Autodock Vina was developed in this section. We predicted whether an active constituent could enter the active pocket of the target protein successfully and calculated the affinities between them. We summarized the affinities of the selected active constituents with interactive residues and the count of hydrogen bonds formed between interactive residues, which are set out in Table3 and Figure 8.Table 3
Molecular interactions of key target and active constituent.
IntegrationTargetsAffinity (kcal/mol)Number of hydrogen bondsH bond interacting residuesQuercetinMYC−7.63GLU957, GLN954, ALA955LuteolinABCG2−82PHE439, THR435QuercetinABCG2−7.81THR435ApigeninABCG2−8.31THR435BaicaleinCDK1−81LYS33QuercetinCDK1−83GLU81, ASP86, GLU12Figure 8
Molecular docking of key target proteins and active constituents. (a) For quercetin and MYC. (b–d) For luteolin, apigenin, and quercetin, respectively. (e) For baicalein and CDK1. (f) For quercetin and CDK1.
(a)(b)(c)(d)(e)(f)


#### 3.6.1. Docking of Quercetin on MYC

As shown in Table3 and Figure 8, the binding affinity of this combination was −7.6 kcal/mol, and GLN958, GLN954, ALA955, MET253, and GLU957 were identified as interactive residues. Quercetin was bound with MYC by forming 4 hydrogen bonds with GLN954, GLN957, and ALA955. The MET253 residue was found to form a pi-sulfur and the ALA955 residue was formed with pi-alkyl. In addition, there were 8 van der Waals interactions between quercetin and ILE961, LYS256, GLU957, TYR252, TYR259, LEU951, GLN954, and GLN958.


#### 3.6.2. Docking of the 3 Active Constituents on ABCG2

As presented in Table3, the binding affinity of the luteolin upon ABCG2 was −8 kcal/mol. The PHE439 residue interacted with luteolin by forming one hydrogen bond and the THR435 residue formed one hydrogen bond. The binding affinity between apigenin and ABCG2 was the highest among 3 integrations. Although only one hydrogen bond was found between apigenin and ABCG2, we observed 8 van der Waals interactions between apigenin and ILE543, MET549, VAL401, PHE432, LEU405, ASN436, and SER440. Interestingly, the PHE439 residue formed 3 pi-pi stacked interactions. As for quercetin, one hydrogen bond, four pi-alkyl interactions, and one pi-anion interactions were identified.


#### 3.6.3. Docking of the 2 Active Constituents on CDK1

The binding affinity of 2 active constituents on CDK1 was −8 kcal/mol. From Figure8, there were 3 hydrogen bonds provided by the GLU81, ASP86, and GLU12 residues in the interaction with quercetin. In addition, the PHE80, PHE82, LEU83, GLY11, GLY13, and ASP146 residues constituted 6 van der Waals interactions. However, there was no hydrogen bond between baicalein and CDK1.


## 3.1. SBH Active Constituent-Target Network

In our study, we obtained 64 active constituents in total. All active constituents were ranked and detailed information of each active constituent is shown in Table1. The result of SBH active constituent-target network is presented in Supplementary Figure S1, which was composed of 441 nodes (including 64 active constituent nodes and 377 target nodes) and 1161 edges, and the result of the top 10 active constituents and targets according to the degree is set out in Supplementary Table S1.Table 1
Active constituents ofScutellariae barbata herba.
RankingActive constituentsMol01(2R)-5,7-Dihydroxy-2-(4-hydroxyphenyl)chroman-4-oneMol0224-Ethylcholest-4-en-3-oneMol035,7,4′-Trihydroxy-6-methoxyflavanoneMol045,7,4′-Trihydroxy-8-methoxyflavanoneMol055-Hydroxy-7,8-dimethoxy-2-(4-methoxyphenyl)chromoneMol066-HydroxynaringeninMol077-Hydroxy-5,8-dimethoxy-2-phenyl-chromoneMol089,19-Cyclolanost-24-en-3-olMol09BaicaleinMol10BaicalinMol11Beta-sitosterolMol12CampesterolMol13CarthamidinMol14Chrysin-5-methyletherMol15CLRMol16Daucostero_qtMol17DinatinMol18EriodictyolMol19LuteolinMol20MoslosooflavoneMol21QuercetinMol22RhamnazinMol23RivularinMol24SalvigeninMol25SitosterolMol26Sitosterol acetateMol27Stigmasta-5,22-dien-3-ol-acetateMol28Stigmastan-3,5,22-trieneMol29StigmasterolMol30WogoninMol31Scutebarbatine FMol326-O-Nicotinoylscutebarbatine GMol33ApigeninMol34Scutelinquanine DMol35Scutebarbatine EMol36Scutebarbatine NMol376-O-(2-Carbonyl-3-methylbutanoyl)scutehenanine AMol386-O-Acetylscutehenanine AMol397-O-Nicotinoylscutebarbatine HMol40Barbatellarine BMol41Barbatin CMol42Barbatin EMol43Scutebarbatine CMol44Scutebarbatine DMol45Scutebarbatine HMol46Scutebata AMol47Scutehenanine AMol48Scutehenanine DMol496,7-Di-O-nicotinoylscutebarbatine GMol506-O-Nicotinoyl-7-O-acetylscutebarbatine GMol51Barbatin AMol52Barbatin BMol53Barbatin DMol54Pheophorbide aMol55Scutebarbatine BMol56Scutehenanine BMol57Scutehenanine CMol58Scutebarbatine OMol59Scutebarbatine GMol60ScutellarinMol617-Acetoxybarbatin CMol626-Acetoxybarbatin CMol63Barbatin FMol64Barbatin H


## 3.2. Potential Active Target Proteins and PPI Network

A total of 377 targets of SBH and 6390 target genes of CRC were obtained, and the intersection of the two groups was regarded as PATPs (Figure2). A total of 297 proteins were identified and used for constructing the PPI network. After limiting the species to Homo sapiens and setting a minimum confidence score as 0.95, 150 proteins were retained as shown in Figure 3. TP53, JUN, AKT1, MAPK1, PIK3CA, etc. were the core proteins of PPI.Figure 2
Potential active target proteins.Figure 3
PPI network of potential active target proteins. The nodes represent proteins, and the color of the node represents the degree value.


## 3.3. GO Enrichment Analysis

We constructed a cluster analysis using MCODE before GO enrichment analysis, and 9 clusters were obtained altogether (Table2). Afterwards, the top 3 clusters with high cluster scores were used for the procedure of GO enrichment analysis as shown in Figure 4. Interestingly, the majority of biological processes in cluster 1 were described as follows: negative regulation of G1/S transition of mitotic cell cycle (45.16%), cyclin-dependent protein serine/threonine kinase regulator activity (25.81%), and regulation of G1/S transition of mitotic cell cycle (19.35%), and they are of great importance to cell proliferation and influence the occurrence of cancer. Cluster 2 was constructed by membrane protein ectodomain proteolysis and notch receptor processing, which is ligand-dependent. Cluster 3 was composed of chemokine-mediated signaling pathway (40.51%), lipopolysaccharide-mediated signaling pathway (16.46%), mononuclear cell migration (12.66%), etc.Table 2
Detailed information of 9 clusters using MCODE.
ClusterScoreNodesEdges17.846145125.661434.211204043.4111753.3334563.33345733382.933162292.8891013Figure 4
GO enrichment analysis of clusters 1–3 (a–c). For the whole figure, different colors represent different biological process groups. The network diagrams represent the connections between biological process groups, the pie charts represent the proportion of each group, and the bar charts represent the detailed term of each group.
(a)(b)(c)


## 3.4. KEGG Pathway Analysis

KEGG pathway analysis was conducted with theClusterProfiler package of R language, and 174 significant pathways (adjusted p value < 0.05) were identified. We sorted them with adjusted p value and the top 30 terms are shown in Figure 5(a). Eight KEGG pathways related to human solid cancers including prostate cancer, small cell lung cancer, pancreatic cancer, bladder cancer, non-small cell lung cancer, hepatocellular carcinoma, colorectal cancer, and breast cancer were gathered. KEGG website summarized 9 pathways associated with CRC, eight of which were significantly enriched in our results (Figure 5(b)). Notably, the result of the colorectal cancer pathway is significant with adjusted p value < 0.05.Figure 5
KEGG pathway analysis of potential active targets. (a) For the top 30 terms sorted by adjustedp value and (b) for the eight pathways related to CRC summarized by KEGG website. The color of each term represents the different adjusted p values <  0.05 and the size of each term represents the gene counts.
(a)(b)


## 3.5. Identification of Key Target Protein and Active Constituent

In this section, we did an analysis of DEGs in CRC firstly using GEO2R. Two datasets including GSE32323 and GSE44076 were screened from the GEO database. A total of 369 DEGs were obtained from the GSE32323 dataset and 541 DEGs were obtained from the GSE44076. We normalized the proteins corresponding to these DEGs by the STRING database. Seventeen key target proteins were remaining after intersecting with PATPs (Figure6): ABCG2, SFRP1, MMP7, MYC, CA2, NR3C2, HSD17B2, PLAU, ADH1C, HSD11B2, MMP3, XDH, FABP6, CDK1, MMP1, MMP12, and SPP1. We constructed a network to show the active constituents that corresponded to these target proteins as presented in Figure 7. Among the active constituents mapped by key target proteins, quercetin, luteolin, baicalein, etc., were identified to be the top 10 active constituents in terms of degree value in the SBH active constituent-target network. After a literature review of these key target proteins, we chose MYC, ABCG2, and CDK1 for molecular docking.Figure 6
Intersections of differentially expressed genes related proteins and potential active target proteins.Figure 7
The network of key target proteins and related active constituents. The circular nodes represent the key target proteins, and the diamond nodes represent the active constituents. The color of the node represents the degree value.


## 3.6. Molecular Docking Analysis

To verify how an active constituent binds to target as previously referred to a molecular docking using Autodock Vina was developed in this section. We predicted whether an active constituent could enter the active pocket of the target protein successfully and calculated the affinities between them. We summarized the affinities of the selected active constituents with interactive residues and the count of hydrogen bonds formed between interactive residues, which are set out in Table3 and Figure 8.Table 3
Molecular interactions of key target and active constituent.
IntegrationTargetsAffinity (kcal/mol)Number of hydrogen bondsH bond interacting residuesQuercetinMYC−7.63GLU957, GLN954, ALA955LuteolinABCG2−82PHE439, THR435QuercetinABCG2−7.81THR435ApigeninABCG2−8.31THR435BaicaleinCDK1−81LYS33QuercetinCDK1−83GLU81, ASP86, GLU12Figure 8
Molecular docking of key target proteins and active constituents. (a) For quercetin and MYC. (b–d) For luteolin, apigenin, and quercetin, respectively. (e) For baicalein and CDK1. (f) For quercetin and CDK1.
(a)(b)(c)(d)(e)(f)


### 3.6.1. Docking of Quercetin on MYC

As shown in Table3 and Figure 8, the binding affinity of this combination was −7.6 kcal/mol, and GLN958, GLN954, ALA955, MET253, and GLU957 were identified as interactive residues. Quercetin was bound with MYC by forming 4 hydrogen bonds with GLN954, GLN957, and ALA955. The MET253 residue was found to form a pi-sulfur and the ALA955 residue was formed with pi-alkyl. In addition, there were 8 van der Waals interactions between quercetin and ILE961, LYS256, GLU957, TYR252, TYR259, LEU951, GLN954, and GLN958.


### 3.6.2. Docking of the 3 Active Constituents on ABCG2

As presented in Table3, the binding affinity of the luteolin upon ABCG2 was −8 kcal/mol. The PHE439 residue interacted with luteolin by forming one hydrogen bond and the THR435 residue formed one hydrogen bond. The binding affinity between apigenin and ABCG2 was the highest among 3 integrations. Although only one hydrogen bond was found between apigenin and ABCG2, we observed 8 van der Waals interactions between apigenin and ILE543, MET549, VAL401, PHE432, LEU405, ASN436, and SER440. Interestingly, the PHE439 residue formed 3 pi-pi stacked interactions. As for quercetin, one hydrogen bond, four pi-alkyl interactions, and one pi-anion interactions were identified.


### 3.6.3. Docking of the 2 Active Constituents on CDK1

The binding affinity of 2 active constituents on CDK1 was −8 kcal/mol. From Figure8, there were 3 hydrogen bonds provided by the GLU81, ASP86, and GLU12 residues in the interaction with quercetin. In addition, the PHE80, PHE82, LEU83, GLY11, GLY13, and ASP146 residues constituted 6 van der Waals interactions. However, there was no hydrogen bond between baicalein and CDK1.


## 3.6.1. Docking of Quercetin on MYC

As shown in Table3 and Figure 8, the binding affinity of this combination was −7.6 kcal/mol, and GLN958, GLN954, ALA955, MET253, and GLU957 were identified as interactive residues. Quercetin was bound with MYC by forming 4 hydrogen bonds with GLN954, GLN957, and ALA955. The MET253 residue was found to form a pi-sulfur and the ALA955 residue was formed with pi-alkyl. In addition, there were 8 van der Waals interactions between quercetin and ILE961, LYS256, GLU957, TYR252, TYR259, LEU951, GLN954, and GLN958.


## 3.6.2. Docking of the 3 Active Constituents on ABCG2

As presented in Table3, the binding affinity of the luteolin upon ABCG2 was −8 kcal/mol. The PHE439 residue interacted with luteolin by forming one hydrogen bond and the THR435 residue formed one hydrogen bond. The binding affinity between apigenin and ABCG2 was the highest among 3 integrations. Although only one hydrogen bond was found between apigenin and ABCG2, we observed 8 van der Waals interactions between apigenin and ILE543, MET549, VAL401, PHE432, LEU405, ASN436, and SER440. Interestingly, the PHE439 residue formed 3 pi-pi stacked interactions. As for quercetin, one hydrogen bond, four pi-alkyl interactions, and one pi-anion interactions were identified.


## 3.6.3. Docking of the 2 Active Constituents on CDK1

The binding affinity of 2 active constituents on CDK1 was −8 kcal/mol. From Figure8, there were 3 hydrogen bonds provided by the GLU81, ASP86, and GLU12 residues in the interaction with quercetin. In addition, the PHE80, PHE82, LEU83, GLY11, GLY13, and ASP146 residues constituted 6 van der Waals interactions. However, there was no hydrogen bond between baicalein and CDK1.


## 4. Discussion

SBH has been used in TCM for thousands of years, the function of which is described as “heat-clearing and detoxifying” in TCM therapeutic principle and theChinese Pharmacopoeia 2015. In spite of that, SBH has been used in various cancers, specially treated with CRC [17].In our study, we investigated the potential mechanisms of SBH on CRC and predicted the pivotal active constituents and target proteins. The SBH active constituent-target network consisted of 64 active constituents and 441 targets, revealing the pharmacological foundation of SBH. The potential active constituents and targets were deduced in this section. The majority of potential active constituents of SBH are flavonoids. Quercetin (degree = 147) is a polyphenolic flavonoid with underlying anticancer activity, which exists ubiquitously in the vegetal food source, especially in various traditional Chinese medicine [51]. Accumulation of in vitro and in vivo studies has concentrated on potential chemopreventive activity and underlying mechanisms of quercetin in CRC. In vitro studies have identified the following effects of quercetin: induction of cell cycle arrest at the G1 phase, induction of apoptosis and autophagy, decreased expression of cyclooxygenase-2 (COX-2), and heat shock protein synthesis (HSP90AA1, degree = 18), etc. [52–57]. In vivo studies have found that quercetin regulated proliferation and apoptosis via suppressed expression of cyclooxygenase-1 (COX-1), COX-2, and inducible nitric oxide synthase (iNOS) [58]. Luteolin (degree = 93) is a naturally occurring flavonoid, which also serves as a dietary flavonoid [59]. The effects of luteolin in CRC have been proved similar to quercetin such as cell growth inhibition and induction of apoptosis [60, 61]. Wogonin (degree = 55) and baicalein (degree = 37) also belong to the flavonoid, and both of them have shown a significant antitumor effect that has been verified experimentally on CRC cells [62]. In addition, apigenin (degree = 40), beta-sitosterol (degree = 37), stigmasterol (degree = 31), etc., which are not flavonoids but are confirmed to have effects on CRC in vivo or in vitro, also contribute to the primary potential active constituents. A total of 441 targets were predicted and we noticed that many targets were contacted with several active constituents such as NR3C1, HSD11B2, PTGS2, CYP23B1, and PTGS1, and we speculated the top 10 targets could be important in the treatment of CRC as shown in Supplementary Table S1. NR3C1, a receptor of glucocorticoids [63], plays an important role in inflammatory responses and cellular proliferation and differentiation [64]. However, NR3C1 (degree = 24) has been identified likely to be a CRC suppressor gene [65]. Moreover, NR3C1 was proved as a differentially expressed gene in CRC by bioinformatics analysis [66]. Prostaglandin G/H synthase 2 (PTGS2, degree = 21) and prostaglandin G/H synthase 1 (PTGS1, degree = 19) are crucial enzymes in the conversion of arachidonate to prostaglandin H2 (PGH2), which are well known as the targets of nonsteroidal anti-inflammatory drugs (NSAIDs) such as aspirin and ibuprofen [67]. Previous research has established that inhibition of the PGHSs with NSAIDs reduced the development of colon cancer [68]. PTGS2-positive patients were faced with an increased risk of CRC recurrence and poorer CRC-specific survival [69]. In particular, PIK3CA, the gene coding for PI3K p110α and the catalytic subunit of PI3K [70], was found in this network. The mutation of PIK3CA exists in approximately 15–20% of CRC [71], which influences the activation of PI3K-AKT signaling pathway and mTOR signaling pathway [72–74]. PI3KCA was simultaneously targeted by 2 active constituents: scutebarbatine E and scutebarbatine N. Scutebarbatine E is a kind of neo-clerodane diterpenoid alkaloid extracted from SBH, and an in vitro experiment has shown significant cytotoxic activities against HT29 CRC cells [75]. Scutebarbatine N belongs to norditerpenoid alkaloids, which also showed significant cytotoxic activities in HT29 CRC cells [76]. However, their specific mechanisms have not been reported yet, so our study contributes to this aspect.We defined the intersection of SBH target proteins and CRC-related proteins as the PATPs. In order to describe the target protein’s function fully, a PPI network was constructed using the PATPs, as could be seen in Figure3. TP53, JUN, AKT1, MAPK1, etc. were the core proteins of PPI. To date, previous studies have revealed that TP53 is a tumor suppressor gene in many tumor types, which induces cell cycle arrest and apoptosis [77]. TP53 mutation is the most typical phenomenon in human cancers [78]. An international collaborative study has reported the occurrence of TP53 mutations in CRC, which was found in 34% of the proximal colon tumors as well as 45% of the distal colon and rectal tumors [79]. Although the mechanisms of TP53 mutation are not understood fully, there is no doubt TP53 has a great capability as a therapeutic strategy in the future. JUN was identified to be bound to the USP28 (a nuclear-localized deubiquitinase related to DNA damage response checkpoint and MYC protooncogene stability) promoter and involved in GTPase KRas (KRAS)-mediated transcriptional activation of USP28 in CRC [80]. AKT1 is one of the AKT kinases, involved in many biological processes such as proliferation, cell survival growth, and angiogenesis, and the therapeutic potential of inhibitors targeting PI3K-AKT pathway in cancer has been discussed [81]. PI3K-AKT pathway plays a pivotal role in the mechanisms of traditional Chinese medicine, which are used frequently when treated with CRC. The majority of them are involved in the inhibition of PI3K-AKT pathway through SRC and AKT1 [82]. Three MAP kinases (MAPK) including MAPK1, MAPK14, and MAPK10 are found in Figure 3, and MAPK plays a crucial role in the MAPK/ERK cascade [83]. Historically, research has investigated that MAPK1 and MAPK14 were associated with cancer risk and survival in CRC [84]. In addition to the 4 proteins mentioned above, other proteins, such as VEGFA, IL6, CDK1, MYC, could also induce the survival of CRC cells. Therefore, we speculated that one of the core functions of the PPI network was involved in the regulation of cancer cells.KEGG pathway analysis was accomplished usingR package ClusterProfiler, and the most significant 30 of 174 pathways are set out in Figure 5. The most striking result emerging from the data is that colorectal cancer pathway (adjusted p value: 4.89E − 15) is significant, and the results also indicate that SBH has the potential to treat diverse cancers such as prostate cancer [85, 86], lung cancer [87, 88], breast cancer [86, 89, 90], and pancreatic cancer [91], which has been confirmed by previous studies. The KEGG website has summarized 9 pathways related to colorectal cancer and we investigated the enrichment information in our study of these pathways. Eight of 9 pathways were found in our study including MAPK signaling pathway (adjusted p value: 1.46E − 12), ErbB signaling pathway (adjusted p value: 2.82E − 11), cell cycle pathway (adjusted p value: 1.18E − 10), p53 signaling pathway (adjusted p value: 7.10E − 17), mTOR signaling pathway (adjusted p value: 0.001478074), PI3K-Akt signaling pathway (adjusted p value: 1.58E − 18), apoptosis pathway (adjusted p value: 1.78E − 16), and Wnt signaling pathway (adjusted p value: 7.62E − 05). It indicates that the therapeutical effects of TCM depend on the cooperation of multiple pathways. We speculated the PI3K-AKT signaling pathway was the most important pathway, which was the most significant among them with a total of 52 genes gathered. With many fundamental cellular functions such as proliferation, growth, and survival identified [92, 93], the potential of PI3K-AKT signaling pathway in cancer has already been discussed [81, 94]. Inhibition of PI3K has become a new therapeutic strategy in CRC with PIK3CA mutation [95]. Interestingly, PI3K-AKT signaling pathway corresponds to mTOR signaling pathway. mTOR, a serine/threonine protein kinase, is known as an important downstream effector of AKT and related to the activation of AKT [96]. In addition, MAPK signaling pathway is interconnected with PI3K-AKT signaling pathway [97], and approximately 30–40% CRC patients harbor a mutation in KRAS [70], which is a part of RAS-RAF-MAPK cascade joining in the cellular function of CRC cells [98]. Therefore, inhibition of both MAPK and PI3K-AKT signaling pathway could be a more effective strategy. Recurrent genetic alterations of Wnt signaling pathway occur in the majority of CRC. The adenomatous polyposis coli (APC) is a tumor suppressor gene and regarded as a central hub in early CRC, and mutations in the APC often dysregulate the Wnt signaling pathway [99]. Our findings indicate that SBH produces the healing efficacy for CRC possibly by regulating the pathways mentioned above.Three key target proteins were selected for molecular docking. MYC is a protooncogene encoding nuclear transcription factor and regulates tumorigenesis [100], which also serves as a downstream effector of Wnt and Ras signaling pathways in CRC [101]. Genetic ablation of c-Myc, a member of the MYC family, suppresses intestinal tumorigenesis in CRC mouse models [102]. Moreover, the inhibitor of c-Myc was identified to be advantageous in anticancer, such as 10058-F4, which could induce apoptosis and differentiation of acute myeloid leukemia cell [103]. Thus, MYC plays an important role in further CRC therapy [104], and in vitro experiment has proved that quercetin induced apoptosis in HT-29 cells and reduced expression of c-Myc [54]. Our molecular docking results indicated that quercetin has a good docking affinity with MYC. The previous set of analyses has identified that quercetin was a pivotal active constituent of SBH, which formed 3 hydrogen bonds with CDK1, and that the binding energy was −8 kcal/mol. Though there is no hydrogen bond found between baicalein and CDK1, baicalein shared the same affinity with quercetin. CDK1 is a serine/threonine kinase which belongs to the CDKs family and regulates the cell cycle, and it has been wildly accepted that CDK1 is the only essential cell cycle CDK [105]. The nonselective inhibitor of CDK-dinaciclib has been shown to arrest cell-cycle progression and inhibit tumor growth [106]. In recent years, literature has identified that CDK1 is a mediator of apoptosis resistance in CRC [107]. Among the three active constituents including quercetin, luteolin, and apigenin targeted ABCG2, apigenin seems to be the most promising active constituent. Multidrug resistance (MDR) is a pivotal factor influencing the efficacy of chemotherapy and the prognosis of tumor patients. Overexpression of adenosine triphosphate (ATP)-binding cassette (ABC) transporter family is one of the most important mechanisms of MDR, and the major ABC transporters include P-glycoprotein (P-gp/ABCB1), breast cancer resistance protein (BCRP/ABCG2), and multidrug resistance-associated protein 2(MRP2/ABCC2) [108]. Importantly, one of the mechanisms of resistance to irinotecan in CRC is that ABCG2 gene encodes ABC efflux transporter and reduces intracellular drug accumulation [109], so targeting ABCG2 is an effective therapeutic principle to enhance the efficacy of irinotecan [110]. Interestingly, it has been reported that decreasing expression of ABCG2 and ABCB5 may induce the depletion of c-Myc and enhance the chemosensitivity of colon cancer stem cells (CSCs) [111]. We did not perform molecular docking analyses for the remaining 14 of the 17 key targets, due to the lack of sufficiently strong evidence for their use as therapeutic targets at this time. The expression level of HSD11B2 gene was significantly increased in CRC tissues, and the ectopic expression of HSD11B2 gene promoted the metastasis of CRC [112]. FABP6, hypermethylated SFRP1, XDH, PLAU, ADH1C, HSD17B2, and SPP1 have been identified more as a colorectal cancer biomarker than as a therapeutic target at present [113–119]. NR3C2, CA2, and MMP1 were identified as key target proteins in another network pharmacological pharmacology analysis of the colorectal cancer. Jin et al. [120] performed molecular docking between these proteins and quercetin, stigmasterol, and baicalein which were also identified in our study as shown in Figure 7. MMP1, MMP3, MMP7, MMP12, and MMP13 belong to the matrix metalloproteinase family, and MMPs were able to participate in the tumor metastasis process by degrading the ECM. However, MMP2 and MMP9 were significantly related to colorectal cancer and could be regulated by Chinese medicine [121]. These targets are closely related to CRC, which may become potentially therapeutic targets and provide reference for future research.


## 5. Conclusions

Traditional Chinese medicine is characterized by multiple components, multiple targets, and multiple pathways. Network pharmacology, along with molecular docking, bioinformation, and system biology, provides an available methodology to uncover the complex therapeutical effects of TCM. Prior studies have confirmed that SBH presents noticeable antitumor effects. We set out to investigate the potential mechanisms of SBH and hunt pivotal active constituents, targets, and pathways.In this study, a total of 64 active constituents of SBH were obtained from 7 TCM databases and literature, and these active constituents were associated with 377 targets, which were mapped to predicted targets of CRC to get 297 common targets treated as PATPs. After that, a PPI network was constructed to demonstrate the interactions between PATPs. The second major investigation was that we did a GO and KEGG analysis using PATPs. Finally, we did a differentially expressed gene analysis of CRC and 229 DEGs were obtained. After the DEGs related proteins were mapped to PATPs, seventeen key target proteins were remaining for molecular docking. The result indicated that quercetin, luteolin, wogonin, and apigenin were the effective active constituents of SBH. NR3C1, PIK3CA, TP53, JUN, MAPK1, MYC, CDK1, and ABCG2 were the major targets of SBH in the treatment of CRC. The most obvious finding emerging from GO analysis was that the core biological process regulated by SBH was the regulation of cell cycle. One of the most significant findings from KEGG analysis was that pathways were significantly enriched in CRC and its related pathways. Molecular docking results reveal that SBH’s active constituents have an acceptable binding affinity with MYC, CDK1, and ABCG2, all of which have shown the potential to treat with CRC. Interestingly, we found that targeting MYC and ABCG2 could contribute to enhancing the efficacy of chemotherapy. A limitation of this study is that further experiments are necessary to demonstrate our findings.



---

*Source: 3905367-2021-08-02.xml*