Based on Network Pharmacology to Explore the Molecular Mechanism of Buzhong Yiqi Decoction for the Treatment of Gastric Cancer

Background Buzhong Yiqi Decoction (BZYQD) has been widely accepted as an alternative treatment for gastric cancer (GC) in China. The present study set out to determine the potential molecular mechanism of BZYQD in the treatment of GC by means of network pharmacology, molecular docking, and molecular dynamics simulation. The potential active ingredients and targets of BZYQD were screened out through the Traditional Chinese Medicine Systems Pharmacology (TCMSP). GC-related targets were screened out through the GeneCards database, and the intersection targets of BZYQD and GC were obtained by using the Venn diagram online tool. Then, the TCM-Active Ingredient-Target network was constructed by using the Cytoscape, and the protein-protein interaction (PPI) network was constructed by using the STRING database. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the effective targets of BZYQD in GC were performed through the Metascape platform. Finally, the molecular docking between the compounds and the target proteins was performed by using the AutoDock Vina software. The simulation of molecular dynamics was conducted for the optimal protein-ligand complex obtained by molecular docking using the Amber18 software.


Introduction
Gastric cancer (GC) represents the most prevalent gastrointestinal malignant tumor worldwide, with its ever-increasing prevalence being very alarming. According to the statistics released in 2021, the morbidity of GC ranks fth and the mortality ranks fourth in the global malignant tumors, causing more than 1 million new cases and 769,000 deaths in 2020 [1] . At present, the pathogenesis of GC has not been thoroughly elucidated. The existing knowledge indicates that bacteria (such as Helicobacter pylori), smoking, and inherited mutations of speci c genes (such as CDH1 gene) may underlie the prominent risk factors for GC [2] . Currently, surgical resection remains the preferred therapeutic method for GC patients, but unfortunately, the majority of GC patients have already progressed into the middle and advanced stages when they are diagnosed, with a 5-year survival rate of less than 30% [3] . Moreover, chemotherapy often elicits gastrointestinal adverse reactions and bone marrow suppression [4] . Against this backdrop, it is of great clinical signi cance to seek alternative or adjuvant treatment methods with few side effects and favorable curative e cacy for GC patients.
Traditional Chinese medicine (TCM) has been practiced for thousands of years in China. Accumulating studies have demonstrated that TCM combined with modern medicine can enhance e cacy and reduce toxicity in the treatment of human malignancies [5] . Buzhong Yiqi Decoction (BZYQD) is derived from the "Spleen and Stomach Theory" written by Li Dongyuan, a famous doctor in the Jin Dynasty. BZYQD is composed of 8 kinds of herbs: Astragalus membranaceus ("Huangqi" in Chinese), Atractylodes macrocephala ("Baizhu" in Chinese), Ginseng ("Renshen" in Chinese), Licorice ("Gancao" in Chinese), Tangerine Peel ("Chenpi" in Chinese), Radix bupleuri ("Chaihu" in Chinese), Angelica sinensis ("Danggui" in Chinese), and Cimicifugae ("Shengma" in Chinese). In recent years, many scholars have unveiled the great implication of BZYQD in the adjuvant treatment of GC, which contributes to reducing side effects, improving quality of life, and prolonging survival. BZYQD can restore impaired immune function, enhance the viability of natural killer cells, and exercise notable anti-tumor effects [6] . Also, BZYQD can inhibit the expression of PD-L1 via the PI3K/AKT pathway, repress the increase of CD8+PD-1+T cells in GC patients induced by chemotherapy, improve immune function, enhance tumor lymphocyte in ltration, and impede the progression of GC [7] . Yu et al. [8] have revealed that BZYQD combined with cisplatin in the treatment of malignant tumors can reverse cisplatin resistance by inducing the accumulation of reactive oxygen species (ROS), and ROS accumulation is known to activate apoptosis and autophagy through oxidative stress. Pei et al. [9] have suggested that BZYQD can improve the immune function of GC patients, facilitate the Th1/Th2 immune imbalance, and reduce the adverse reactions of chemotherapy.
Network pharmacology, an emerging discipline based on system biology and computer technology, has been extensively applied in the research eld of TCM compounds in recent years due to its functionality of predicting effective ingredients, targets, toxic, and side effects of drugs [10,11] . Therefore, this study set out to screen out the active components of BZYQD in the treatment GC and predict the related targets and signaling pathways with the support of the network pharmacology to explore the mechanism of action of BZYQD in the treatment of GC. The detailed owchart of the present study is shown in Figure 1 In line with the selection criteria of oral bioavailability (OB) ≥ 30% and drug-likeness (DL) ≥ 0.18, the active ingredients of 8 herbs in BZYQD were searched in the Traditional Chinese Medicine Systems Pharmacology (TCMSP) [12] (https://tcmspw.com/tcmsp.php), and then the targets corresponding to each active ingredient was searched. Finally, the collected target names were converted into the ID names of the corresponding genes by using the Uniport protein database [13] .
2.2 Prediction of the potential targets of BZYQD in the treatment of GC.
With "gastric cancer" as the keyword and Score > 10 as the screening criteria, GC-related targets were retrieved from the GeneCards database [14] . The targets of BZYQD and the targets of GC were analyzed by using the Venn diagram online tool, and a Venn diagram was made. The intersection between the two groups of targets was considered the potential targets of BZYQD against GC.
2.3 Construction of a "TCM-Active Ingredient-Target" network.
The potential targets of BZYQD in the treatment of GC and the corresponding active ingredients were imported into Cytoscape 3.8.0 [15] to construct a "TCM-Active Ingredient-Target" network to clarify the pharmacological mechanism of BZYQD against GC. The degree value re ected the importance of the node in the network. The higher the degree value, the more important the node. The key active ingredients were screened out according to the degree value.

Construction of protein-protein interaction (PPI) network and cluster analysis.
The target proteins corresponding to the components of BZYQD in GC were imported into the STRING [16] (https://string-db.org) to perform PPI analysis, with the species set to "Homo sapiens", the con dence set to > 0.9, and unconnected nodes hidden. Then, the tsv le of PPI analysis results was imported into Cytoscape to construct a PPI network. Some target proteins had the same or similar characteristics and functions in cell biological activities, and these target proteins were regarded as a cluster. Proteins in the same cluster are generally thought to play a synergistic role in disease progression [17] . Cluster analysis of target proteins in the bioinformatics network was performed by using the MECODE plug-in of Cytoscape.

Enrichment analysis.
The Metascape platform [18] was used to perform Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on the intersection targets of BZYQD and GC. The results of enrichment analysis were saved and nally visualized by using R software. The biological process (BP), cell component (CC), and molecular function (MF) results of GO enrichment analysis were drawn into bar graphs, and the KEGG pathway enrichment analysis results were drawn into a bubble chart for the further exploration of the biological processes and metabolic pathways of BZYQD in GC.

Molecular docking.
The top ve components in degree value in the "TCM-Active Ingredient-Target" network were molecularly docked with the top ve core target proteins in degree value in the PPI network. The structure of the compounds was obtained from the PubChem database [19] and the structure of the proteins was from the RCSB PDB database [20] . Then, the AutoDock software was used for molecular docking. Three conformations with the best docking binding energy were selected for the docking binding mode analysis. Discovery Studio Visualizer 2019 software was used for force analysis and PyMOL software was used for mapping.

Molecular dynamics simulation.
The amber18 software [21] was used to perform molecular dynamics simulation on the protein(target)and small molecule ligand (compound) obtained by molecular docking. FF144SB and general Amber force eld (GAFF) parameters were used for proteins and small molecule ligands, respectively, and the AM1-BCC atomic charge was calculated using the Antechamber module. The protein-ligand complex was loaded into the leap module, and hydrogen atoms and antagonist ions were automatically added to neutralize the charge. The TIP3P explicit water model was selected and periodic boundary conditions were set. The work ow of molecular dynamics simulation included energy minimization, heating, equilibration and production dynamics. Heavy atoms of proteins (small molecules) were constrained, and water molecules were subjected to 10,000 steps of energy minimization (including 5000-step steepest descent method and 5000-step conjugate gradient method). The system was slowly heated to 300 K within 50 ps and equilibrated for 50 ps under the NPT ensemble. Molecular dynamics simulation was conducted on the system under the NPT ensemble for 50 ns (25,000 steps in total) with a time step of 2 fs. The trajectory data were saved every 10 ps, followed by analysis using the CPPTRAJ module. The free energy of binding for ligands and proteins was calculated using the MMPBA.py module.

Information on Potential Active Ingredients and Targets in BZYQD.
With OB ≥ 30% and DL ≥ 0. 18  3.2 Information on GC target genes and acquisition of "TCM-Disease" intersection targets.
A total of 12136 GC targets were obtained from the GeneCards database. With a score >10 as the screening criterion, a total of 1383 GC targets were retained. Finally, 136 TCM-Disease intersection targets were obtained using Venn diagram online tool (Figure 2).
The targets and active ingredients of BZYQD in GC were imported into Cytoscape to construct a "TCM-Active Ingredient-Target" network ( Figure 3). It was found that the top vecompounds with the highest degree value in BZYQD were quercetin, kaempferol, nobiletin, naringenin, and formononetin, indicating that these vecompounds might be the effective components of BZYQD in the treatment of GC (Supplementary Table S1).

Construction of PPI Network and Cluster analysis.
The intersection targets were imported into the STRING database for PPI analysis, and the interaction threshold was set to "0.9" and the unconnected nodes were hidden. Then, the PPI analysis results were imported into Cytoscape software to draw a PPI network diagram ( Figure 4). The PPI results demonstrated that the top ve core target genes in degree value were AKT1 (RAC-alpha serine/threonineprotein kinase, Degree = 45), STAT3 (Signal transducer and

GO enrichment and KEGG pathway enrichment analysis.
GO and KEGG enrichment analysis on the targets of BZYQD in GC were performed through the Metascape platform. The results showed that there were 2208 BP, 77 CC, and 145 MF in the GO enrichment analysis. KEGG pathway enrichment analysis yielded 172 signaling pathways. According to the P-value, the top 10 BP, CC, and MF were selected for visualization ( Figure 6). The top 20 pathways of KEGG enrichment analysis were visualized using R software to draw a bubble diagram (Figure 7). The apoptosis signaling pathway, PI3K/Akt signaling pathway, proteoglycan in cancer, IL-17 signaling pathway, TNF signaling pathway, and HIF-1 signaling pathway were mainly involved (Supplementary  3.6 Molecular docking results. The ve compounds (quercetin, kaempferol, nobiletin, naringenin, and formononetin) with the highest degree value were docked with ve core target proteins (AKT1, STAT3, TP53, MAPK1, and MAPK3) in the PPI network by using AutoDock software (Supplementary Table S3). The results demonstrated that the key compounds showed good binding activity with the core target proteins, and the best binding energy was MAPK3-naringenin (-8.08 kcal/mol). The molecule formed two hydrogen bonds with the protein amino acid residues Ser74 and Tyr81 at a distance of 2.57 Å. Moreover, it formed a hydrophobic interaction with Gly54, Ile73, Pro75, and Tyr81, followed by MAKP3-formononetin (-7.97 kcal/mol) and MAPK1-kaempferol (-7.92 kcal/mol). The molecular docking diagram is shown in Figure 8.

Molecular dynamics simulation results.
Based on the binding energy results, molecular dynamics simulation was performed for MAPK3-Naringenin. The Root-mean-square deviation (RMSD), Radius of gyration (Rog), Root-mean-square uctuation (RMSF) curve, and binding free energy were calculated.
The RMSD curve represents the variations in protein conformation. The results showed that the RMSD uctuated greatly in the rst 30ns and then gradually stabilized, indicating that the protein-ligand complex formed a stable binding conformation.
The Rog curve displays the compactness of the overall protein structure. It could be seen from the results that the radius of gyration was stable, indicating the compact structure of the protein during the kinetic process and the unchanged conformation.
The RMSF curve represents the variations in the conformation of amino acid residues. It could be seen from the results that the C-terminal region of the protein had greater residue exibility than other regions.
Since this region was a loop-based polypeptide chain, it was more exible. The results of RMSD, Rog, and RMSF are shown in Figure 9.
The trajectory of the RMSD plateau (30-50 ns) was calculated, allowing the calculation of the free energy of binding for ligand-protein interactions based on the MMGBSA equation. The results showed that the binding free energy of MAPK3-Naringenin was -25 kcal/mol (Supplementary Table S4). Van der Waals potential (-32 kcal/mol) and electrostatic interaction (-26 kcal/mol) were conducive to the combination of the two, while polar solvation (EGB, 37 kcal/mol) was not conducive to the interaction of the two.
The binding site of MAPK3-Naringenin is a binding site with mixed hydrophobic and hydrophilic properties. Residues forming hydrophobic accumulation with ligands included Leu150, Leu101, etc., and residues forming electrostatic interaction included Lys48, Asp161, etc. Hydrophobic accumulation and electrostatic interaction (including hydrogen bonding) were crucial for stable ligand binding. The MAPK3-Naringenin interaction diagram is shown in Figure 10.

Discussion
Through the network pharmacological method, this study elucidated that quercetin, kaempferol, nobiletin, naringenin, and formononetin were the active ingredients of BZYQD in the treatment of GC. Prior studies have shown that quercetin can regulate the PI3K/Akt pathway by affecting AKT phosphorylation, downregulating AKT kinase, and inhibiting EMT in GC cells, thereby repressing GC cell migration and invasion [22][23] . Kaempferol can promote autophagy of GC cells, increase the conversion from LC3-I to LC3-II, and down-regulate the expression of p62 protein in GC [24] . Nobiletin triggers endoplasmic reticulum stress-mediated apoptosis and autophagy of GC cells via the downregulation of the Akt/mTOR signaling [25] . Naringenin can inhibit the proliferation, migration, and invasion of GC cells, and induce apoptosis by inhibiting the phosphorylation of AKT and decreasing the expressions of its downstream target molecules [26] . Formononetin can depress the growth and migration of GC cells by regulating the Wnt/β-Catenin and AKT/mTOR pathways [27] .
The results of PPI analysis showed that the targets of BZYQD relating to the treatment of GC mainly involved AKT1, STAT3, TP53, MAPK1, and MAPK3. Akt is a protein kinase that regulates the cell cycle and apoptosis. Activated Akt can interact with several key downstream effectors and phosphorylate, thereby suppressing cell apoptosis [28] . STAT3, a transcription factor of the JAK-STAT signaling pathway, has the ability to promote GC cell proliferation, invasion, and migration [29] . In the tumor microenvironment, tumor cells can secrete interleukin, which in turn activates the STAT3 signaling pathway, leading to the progression and deterioration of GC cells [30] . TP53 is commonly accepted as a tumor suppressor gene. The expression of TP53 is increased continuously during the deterioration of normal gastric mucosa into intestinal metaplasia and then into GC, and TP53 mutation can be frequently observed in GC patients [31] .
MAPK is widely implicated in numerous vital biological processes such as cell proliferation, differentiation, and apoptosis through phosphorylation of nuclear transcription factors and related enzymes [32] . MAPK1 activation is indicative of the invasiveness of GC. Down-regulation of MAPK1 expression can retard GC cell growth, promote apoptosis, and suppress migration and invasion [33] . Kim et al. [34] have revealed that MAPK3 expression can be regarded as an independent prognostic index for GC patients receiving resection treatment.
GO enrichment analysis indicated that BZYQD treated GC through a variety of pivotal biological processes such as oxidative stress, chemical stress, lipopolysaccharide reaction, and apoptosis. KEGG pathway enrichment analysis exhibited that the main pathways of BZYQD in the treatment of GC included the apoptosis signaling pathway, PI3K/Akt signaling pathway, proteoglycan in cancer, IL-17 signaling pathway, TNF signaling pathway, and HIF-1 signaling pathway. Among these signaling pathways, apoptosis is undoubtedly highlighted as the most critical signaling pathway in malignant tumors, so the induction of tumor cell apoptosis has emerged as a research hotspot in the treatment of GC [35] . The PI3K/Akt signaling pathway is proposed to be a carcinogenic signaling pathway, which participates in the modulation of diverse biological processes of cancer cells, such as proliferation, apoptosis, chemotherapy resistance, and metastasis [36] . Lee et al. [37] have clari ed that inhibition of the PI3K/Akt signaling pathway can lead to G2/M phase arrest, autophagy, and apoptosis of GC cells.
Proteoglycan, an important part of extracellular matrix proteins, can mediate the migration and proliferation of tumor cells [38] . The IL-17 signaling pathway can dramatically alter the tumor microenvironment by enhancing the secretion of chemokines and cytokines, thereby promoting tumor occurrence and development [39] . The TNF signaling pathway regulates in ammatory responses in the body. Excessive TNF in ammatory cytokines in the tumor microenvironment can elicit tumor growth and destroy the innate immune response of tumor cells [40] . The HIF-1α signaling is regarded as an initiator of GC progression, which can promote GC cell proliferation and angiogenesis, induce EMT, and enhance the invasion and metastasis of GC, and consequently, inhibition of the HIF-1α signaling can be a promising strategy to prevent GC cells from in ltration and metastasis [41] .

Conclusion
Through network pharmacology, we revealed the signi cant advantages of BZYQD in the treatment of GC. BZYQD mainly acts on multiple target genes such as AKT1, STAT3, TP53, MAPK1, and MAPK3 through a variety of compounds including quercetin, kaempferol, nobiletin, naringenin, and formononetin, then interferes with the apoptosis signaling pathway, PI3K/Akt signaling pathway, proteoglycan in cancer, IL-17 signaling pathway, TNF signaling pathway, and HIF-1 signaling pathway, and eventually exerts a therapeutic effect on GC. Figure 1 The detailed owchart of the present study.  Venn diagram of intersection between BZYQD (blue) and GC (red) targets.  The protein-protein interaction network. The sizes and colors of the nodes and edges are illustrated from big to small and blue to orange in descending order of degree values.  KEGG pathway enrichment analysis of the targets of BZYQD in GC. Each bubble represents a KEGG pathway on the vertical axis. The size of the bubble represents the number of genes enriched in each KEGG pathway, and the larger the bubble, the more genes in the pathway. The color of each bubble is determined by the P-value, and the redder the color of the bubble, the smaller the P-value.