Dynamic analysis of a predator-prey state-dependent impulsive model with fear effect in which action threshold depending on the prey density and its changing rate

1.   Introduction

The Oddi sphincter (SO) is a precision smooth muscle device at the junction of the bile duct, pancreas and duodenum. It forms a temporal transitional compound movement (MMC) under the multiple and complex regulations of nerves, humors and local reflexes, which play important roles in controlling bile and pancreatic fluid discharge and preventing reflux ^[1]. Sphincter of Oddi Dysfunction (SOD) is a group of diseases caused by abnormal diastolic function of SO. It is clinically divided into two types, biliary type and pancreatic type, and the former is common ^[2]. SOD often has no evidence of organic change, but it will bring long-lasting and unbearable trouble to patients. In severe cases, it can be secondary to liver damage, abnormal trypsin and even acute pancreatitis, seriously endangering the life and health of patients ^[3]. Therefore, it is of great significance to explore the pathological mechanism and drug treatment of SOD.

Shao Yao Gan Cao Tang (SYGC), sourced from Shang Han Za Bing Lun in 210 CE, is made up of 2 traditional Chinese medicines: Paeoniae radix and Glycyrrhizae radix (1:1), which has been used to treat general muscle pain or tremor in skeletal muscles ^[4]. The Paeoniae Radix can nourish blood and relieve the depressed liver, and Glycyrrhizae Radix can strengthen the spleen and Qi ^[5]. It is reported that SYGC can reduce abdominal pain and muscular cramps ^[4] and suppress duodenal peristalsis during endoscopic retrograde cholangiopancreatography (ERCP) ^[6]. Our clinical study has indicated that SYGC had good curative effects on SOD, relieving abdominal pain symptoms, improving liver function and reducing bile excretion time ^[7]. However, its chemical profile responsible for the therapeutic effects on SOD is still unclear.

Recently, network pharmacology and transcriptomics have been widely used to explore the active components and potential mechanisms of traditional Chinese medicine (TCM) ^[8,9]. In this study, to enable a full assessment of transcriptomics changes in SOD and to increase our understanding of the SYGC mechanism on SOD, ultra-high-performance liquid chromatography coupled with Quadrupole Exactive-Orbitrap high-resolution mass spectrometry (UHPLC-Q Exactive-Orbitrap HR-MS) was first used to identify chemicals of SYGC. Then, transcriptomics and network pharmacology were applied to uncover the mechanism of SYGC in the treatment of SOD. Finally, the binding affinities between the active compounds and the key targets were determined via molecular docking. Figure 1 shows the flowchart of the study design. This study provides a basis for the clinical application of SYGC in the treatment of SOD.

2.   Materials and methods

2.1.   Ultrahigh-performance liquid chromatography quadrupole-Orbitrap high-resolution mass spectrometry (UHPLC-Q-Orbitrap-HR-MS) analysis and identification of target genes of SYGC

The ingredients of SYGC samples were measured by UHPLC-Q-Exactive Orbitrap analysis ^[10] which was performed by an UHPLC system (UltiMate 3000 RS, Thermo Fisher Scientific, USA) coupled with a Q Exactive Orbitrap (QE, Thermo Fisher Scientific, USA) equipped with an electrospray ionization source (ESI). The protonated molecular weights of all identified compounds were calculated within an error of 10 ppm. Following careful comparisons with the retention times and MS/MS spectra of the reference standards, reference literature, Chemical Book and self-built databases, a total of 111 chemicals were identified or tentatively characterized from SYGC. Then, active compounds were screened based on the parameters including Lipinski's "rule of five, " Ghose #violations and GI absorption, which was performed by a SwissADME (http://www.swissadme.ch/) ^[11]. The target information of each component of SYGC was obtained from HERB (http://herb.ac.cn/) ^[12], SwisstargetPrediction (http://swisstargetprediction.ch/) ^[13] and BATMAN-TCM (http://bionet.ncpsb.org.cn/batman-tcm/index.php/Home/Index/index) ^[14].

2.2.   Establishment of SOD model

The animal experiments were performed by the experimental animal welfare ethics review committee of the Shanghai University of TCM. Forty guinea pigs were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China; qualified production number P2020-052). The guinea pigs were fed under a 12 h cycle of light/dark in IVC conditions and had free access to food and water. The guinea pigs were divided into four groups: control group (N), SOD group (M), SYGC gavage treatment group (G), IRE1 inhibit treatment group (IR) as a positive control group (n = 10/group). The M, G and IR groups were injected intravenously with morphine injection at 0.6 mg/kg body weight, and the normal group was injected intravenously with equal volume of normal saline three times a week for a total of 4 weeks. G group was given orally at 12.5 g/kg per day. The IR group was administered by injection of 30 mg/kg twice a week. After the last administration, all animals were fasted for 12 hours, and the abdominal cavity was opened after anesthesia. After the duodenum was dissected, the white papillary protrusion and the texture of the Oddi sphincter tissue were isolated for follow-up experiments. Blood samples were collected for measurement of alanine aminotransferase (ALT) and aspartate aminotransferase (AST) activities.

2.3.   Histological staining

Oddi sphincter tissue samples were fixed in 10% neutral formalin for 48 hours, then dehydrated by conventional gradient, embedded in paraffin, sliced, baked at 60 ℃ for 40 minutes, stained with hematoxylin-eosin (HE) and sealed. The sections were stained with HE to observe the morphological changes such as inflammatory cell infiltration and sphincter injury of Oddi sphincter.

2.4.   Total RNA extractions and RNA-Seq analysis

Total RNA in sphincter tissues from the SOD (n = 3) and control samples (n = 3) were extracted using TRIzol Reagent (TIANGEN, China). The purity, concentration and integrity of the total RNA samples were checked for further analysis, and Samples with RNA integrity number (RIN) ≥ 7 were considered to be of high quality. A transcriptome sequencing using the Illumina sequencing platform (HiSeqTM 2500) was conducted on each total RNA sample by OE Biotech Co., Ltd. (Shanghai, China). The raw data were shown in Supplementary data. The differentially expressed genes (DEGs) were screened with |log₂ fold-change (FC)| ≥ 1.0 and q ≤ 0.05, which applied for gene ontology (GO) and KEGG enrichment analysis.

2.5.   Identification of genes targeted by SYGC in SOD and functional enrichment

To obtain the intersection targets, the key targets of SYGC and SOD DEGs was plotted by https://www.bioinformatics.com.cn, a free online platform for data analysis and visualization. Metascape platform (https://metascape.org/gp/index.html) ^[15] was applied to perform Gene ontology (GO) and pathway enrichment analysis (Wiki, Reactome and KEGG pathway) of the above intersection targets.

2.6.   Construction of a protein-protein interaction network of the intersecting target genes

The STRING (https://string-db.org/) database was applied to create a protein-protein interaction (PPI) network of intersecting target genes. The Cytoscape 3.9.0 (https://cytoscape.org/) was used to visualize complex relationships between the active chemical components and the target genes. Next, Molecular Complex Detection (MCODE) ^[16] and ClusterOne analysis ^[17] were used to find the core targets. MCODE scores ≥ 3 and (node ≥ 3 and P < 0.05) were set as the criteria. The number of nodes ≥ 3 and P < 0.05 were set as the criteria for ClusterOne analysis.

2.7.   Molecular docking validation

The Protein Data Bank (PDB, https://www.rcsb.org/) was used to obtain the crystal structures of core targets. The three-dimensional structures of active ingredients were downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/). Molecular docking was performed to calculate the binding affinity between active ingredients and core targets by the AutoDock Vina (http://vina.scripps.edu/). The corresponding PDB codes were 4af3, 6hky and 4cik for AURKB, KIF11 and PLG, respectively. We first removed the proteins' water molecules, added polar hydrogen. Then, active pockets were built according to the position of ligand in the PDB complex. The binding energy of ligands in PDB structures was applied for positive control.

2.8.   Verifying of hub genes by RT-qPCR

Total RNA was extracted from each sphincter of Oddi tissue using TRIzol (Invitrogen Corporation, CA, USA). cDNA was prepared through cDNA Synthesis SuperMix. Using Gapdh as an internal reference, RT-qPCR was performed to detect the mRNA expression levels of intersecting genes among MCODE genes, ClusterONE genes and genes of significant enrichment pathways. The amplification primers were synthesized by Shanghai Sangon Biological Engineering Technology, as shown in Table 1. The mRNA expression level was normalized to that of Gapdh in the same sample. The relative expression of each target gene was calculated by the 2^-ΔΔCt method.

3.   Results

3.1.   The effect of SYGC in guinea pig with SOD

After the model was successfully constructed, the SYGC and positive drug (IRE1 inhibitor) were administered for four weeks. The control group was given the same dose of physiological saline solution. Compared with the normal group, the sphincter tissue of the model group showed edema, increased inflammatory cell infiltration (Figure 2) and decreased c-kit expression (Figure 3), and the ring muscle had partial disorders and irregularities. Also, symptoms such as submucosal fibrous tissue hyperplasia and smooth muscle hypertrophy appeared. After treatment with SYGC or IRE1 inhibitor, the sphincter tissue showed increased c-kit expression and decreased inflammation infiltration, and ring muscle disorders were reduced, suggesting an improved Cajal cell activity in the Oddi sphincter. In addition, compared with the normal group, the serum ALT and AST levels in the model group were significantly higher (P < 0.01). Compared with the model group, the serum ALT levels in G and IR group were significantly decreased (P < 0.01) (Figure 3), indicating that Shaoyao Gancao Decoction can improve the liver function of SOD guinea pigs.

3.2.   Identification of compounds and related pharmacological parameters of SYGC

To identify the major chemical components, the SYGC samples were analyzed using the UHPLC-Q-Exactive Orbitrap HR-MS analysis. As shown in Figure 4 and Table 2, 23 compounds were identified under in positive ion mode and 88 compounds were identified under in negative ion mode. Then, Swissadme was applied to explore the pharmacological parameters of these compounds. Finally, 32 candidate compounds passed the parameters (Lipinski's "rule of five, " Ghose #violations and GI absorption) and were selected for further research (Table 3). Specifically, 6 ingredients were from Paeonia lactiflora Pall., 22 from Glycyrrhiza uralensis Fisch., and 4 from both herbs.

3.3.   Identification of target genes of SYGC

The target information of 32 active ingredients was obtained from HERB, SwisstargetPrediction and BATMANTCM. Once duplicate genes were deleted, a total of 1023 targets for SYGC were obtained.

3.4.   Identification of SOD differential expression genes and functional enrichment analysis

Furthermore, we explored the dysfunctional genes and pathways in guinea pig SOD using RNA-Seq of sphincter tissues from the Control and SOD groups. The RNA from the three replicate samples from the control and SOD groups was sequenced. In all, 16,281 genes were identified (Supplementary data). To determine the differentially expressed genes (DEGs), a q-value < 0.05 was used as the cut-off value for gene expression in the control and SOD groups using DESeq2. As a result, 649 DEGs including 247 up-regulated and 402 down-regulated genes were screened (Figure 5A). The top 10 enriched GO terms in each category as cellular component (CC), molecular function (MF), biological process (BP) and of the identified DEGs are shown in Figure 5B. The GO terms showed that DEGs were mainly related to negative regulation of cell fate commitment, extracellular space, serine-type endopeptidase inhibitor activity, etc. On the other hand, the results of KEGG enrichment analysis showed that the top enriched KEGG terms were, for example, complement and coagulation cascades, B cell receptor signaling pathway, primary immunodeficiency NF-kappa B signaling pathway (Figure 5C).

3.5.   Determination of SOD-related genes targeted by SYGC and functional enrichment analysis

To identify the intersecting genes between SYGC and SOD, a Venn analysis was performed on the target genes of SYGC and SOD DEGs. As shown in Figure 6A, 52 genes were identified as intersecting genes of SYGC and SOD. Then, these intersecting genes were imported into the Metascape database to carry out GO enrichment analysis and pathway enrichment analysis (Figure 6B). BP terms were mainly found in regulation of ion transport, response to xenobiotic stimulus and leukocyte migration. CC terms were mainly enriched in ion channel complex, side of membrane, and perinuclear region of cytoplasm. MF terms were mainly present in kinase binding, kinase binding-membrane spanning protein tyrosine kinase activity and carbonate dehydratase activity. These factors can exert therapeutic effects on SOD.

In addition, the results of the pathway enrichment analysis mainly involved the B cell receptor signaling pathway, complement system, signaling by receptor tyrosine kinases, Interleukin-4 and Interleukin-13 signaling, as well as muscle contraction, among others (Figure 6C).

3.6.   Analysis of protein-protein interaction network of gene intersection and construction of a regulatory network of targets for the treatment of SOD with SYGC

The STRING database (http://www.string-db.org) was used to investigate the target genes' interactions. There were 57 nodes and 158 edges in the protein-protein interaction network. Then, the complex interactions between active components and potential target genes were visualized using the Cytoscape, including 82 nodes and 214 edges (Figure 7A). Four significant clusters were obtained from ClusterONE analysis (node ≥ 3 and P < 0.05). Cluster 1 consisted of 10 nodes (Figure 7B); Cluster 2 consisted of 11 nodes (Figure 7C); Cluster 3 consisted of 8 nodes (Figure 7D); Cluster 4 consisted of 3 nodes (Figure 7E). Four significant modules were obtained from MCODE analysis (score ≥ 3). Module 1 (score: 8.75) consisted of 9 nodes (Figure 7F); Module 2 (score: 7.125) consisted of 17 nodes (Figure 7G); module 3 (score: 3.6) consisted of 6 nodes (Figure 7H); module 4 (score: 3) consisted of 3 nodes (Figure 7I). Finally, 20 intersecting genes between MCODE genes and ClusterONE genes were obtained, including SERPINE1, MMP9, PLG, CCNB1, CACNA2D2, RAD51, CHRNB4, AURKB, TOP2A, TYMS, KIF11, KCNQ1, CASQ2, BARD1, CHRNA4, IGFBP1, TNNT2, SCN4A, MKI67, CMA1.

3.7.   Molecular docking validation

The present study examined the intersecting genes between MCODE and ClusterONE genes involved in the B cell receptor signaling pathway (AURKB, KIF11) and complement system (PLG). Each of the enriched components is docked with the three genes. The binding energy of ligands in PDB structures was used for positive control (VX6 for AURKB: -8.8 kcal/mol, GCE for KIF11: -9.5 kcal/mol, XO3 for PLG: -7 kcal/mol). Glycycoumarin and licoflavonol exert lower score than VX4 ligand for AURKB target, indicating a strong binding activity. Rhamnocitrin shows a relative lower score compared to GCE ligand for KIF11 target, indicating a good binding activity. Echinatin, homobutein and licoflavonol present a relative lower score compared to XO3 ligand for PLG target, suggesting a good binding activity (Table 4).

Specifically, AURKB showed 12 interactions with glycycoumarin, including unfavorable donor- donor, Pi-sigma, Pi-alkyl and Alkyl, which were connected with GLU 161, LEU 83, ALA 157, VAL 91 and PHE 88 (Figure 8A). It showed 13 interactions with licoflavonol including carbon hydrogen bond, conventional hydrogen bond, unfavorable donor- donor, Pi-sigma, Pi-alkyl and Alkyl, which were connected with VAL 91, LYS 106, ALA 157, GLU 161 and GLY 160 (Figure 8B). KIF11 showed 8 interactions with rhamnocitrin, including conventional hydrogen bonds, Pi-sigma, Pi-alkyl and Pi-anion, which were connected with TRP 127, PRO 137, ALA 133, GLU 116 and ARG 119 (Figure 8C). PLG showed 6 interactions with echinatin including carbon hydrogen bond, conventional hydrogen bond, unfavorable donor- donor, Pi-cation, Pi-Pi stacked and Pi-Pi T-shaped, which were connected with TYR 72, ASP 55, ARG 71, TRP 62 and ARG 35 (Figure 8D). It showed 6 interactions with homobutein including carbon hydrogen bond, Pi-donor hydrogen bond, Pi-Pi stacked, Pi-Pi T-shaped, Pi-anion, Pi-alkyl and Pi-cation, which were connected with ARG 35, TRP 62, TYR 72, ASP 55 and ARG 71 (Figure 8E). It showed 9 interactions with licoflavonol including carbon hydrogen bond, Pi-Pi stacked, Pi-Pi T-shaped, Pi-anion, Alkyl and Pi-cation, which were connected with ASP 55, TRP 62, TYR 72, ASP 57 and ARG 35 (Figure 8F).

3.8.   Verifying of hub genes

To verify the results of bioinformatics analysis, we obtained the top eight genes (Prkcb, Alox5, Adora2a, Lck, Jak3, Rasgrp3, Adipoq and Gh1) involved in the 52 gene-related signal pathways and 14 intersecting genes (Serpine1, Mmp9, Plg, Ccnb1, Cacna2d2, Rad51, Chrnb4, Aurkb, Top2a, Tyms, Kif11, Kcnq1, Casq2 and Bard1) among MCODE genes, ClusterONE genes and genes of significant enrichment pathways. We detected the mRNA expression levels of these genes by RT-qPCR in the Oddi sphincter tissues. Compared with the control group, the M group showed decreased expression of Ccnb1, Mmp9, Rad51, Top2a, Tyms, Kif11, Rasgrp3, Prkcb, Lck, Jak3, Adora2a, Aurkb and Gh1 and increased expression of Cacna2d2, Chrnb4, Alox5 and Plg. Moreover, the expression of Ccnb1, Cacna2d2 and Chrnb4 returned to normal after SYGC treatment (Figure 9).

4.   Discussion

SOD is a key secondary pathological change in the context of gallbladder- and pancreas-related inflammatory diseases and seriously impacts patients' quality of life ^[18,19]. At present, the research on its pathogenesis and effective treatment are in the preliminary stages. TCM has been widely applied for the discovery of candidate drugs ^[20]. SYGC is used for the treatment of pain-related diseases with reducing muscle tension, relieving spasms and providing analgesia ^[21]. Moreover, paeoniflorin, an extract of Shaoyao, can relax the SO muscle via reducing calcium ion influx ^[22]. Isoliquiritigenin, a flavonoid from licorice, relaxed guinea-pig tracheal smooth muscle through the cGMP/PKG pathway ^[23]. Consistent with our clinical study ^[7], we also found that SYGC administration can repair the structure and ultrastructure of the SO in vivo with decreased inflammation infiltration and ring muscle disorders. However, the detailed regulatory mechanism of SYGC action against SOD requires further investigation.

Therefore, we used the systematic pharmacological method to discover the potential molecular mechanisms of SYGC on SOD. At first, a total of 649 DEGs were identified in SOD, which were mainly enriched in complement and coagulation cascades, B cell receptor signaling pathway, primary immunodeficiency and NF-kappa B signaling pathway. Recently, it was well established that immune disorders and inflammation played an important role in the occurrence and development of numerous diseases, including SOD ^[24,25,26]. The B cell receptor (BCR) signaling pathway was crucial for normal B cell development and adaptive immunity ^[27], and B cell-derived IgE may lead to smooth muscle contraction induced by the degranulation of mast cells ^[28]. NF-κB showed a key role in various biological processes, including inflammation, immune response, cell growth and survival and development ^[29,30,31]. NF-κB signaling impeded the recovery of skeletal muscle function after damage ^[32]. In addition, NF-kappaB activation served as a survival factor in B cell, which prevented cell apoptosis ^[33,34]. These results suggest that SOD may exert a dysfunction crosstalk between B cell receptor and NF-κB signaling pathway.

Then, does SYGC ameliorate SOD by regulating the B cell receptor signaling pathway? We first identified 32 candidate compounds using UHPLC-Q-Orbitrap-HRMS and network pharmacology analysis. Interestingly, these chemicals had an effect in the B cell receptor signaling pathway. Additionally, SYGC may improve SOD through multiple pathways including the complement system, Interleukin-4 and Interleukin-13 signaling and muscle contraction. Therefore, these results suggest that the B cell receptor signaling pathway may play a central role in these multiple pathways. For example, IL-4 and IL-13 exerted their signaling action by IL-4Rα/IL-13Rα complexes ^[35], and IL-4 was demonstrated to regulate B-cell receptor signaling in chronic lymphocytic leukemia ^[36]. The complement system, an essential contributor of innate immunity, was also important in regulating B cell responses at multiple stages of the peripheral response ^[37]. Combined with the effect of SYGC on muscle cramps ^[4,38], we speculate that SYGC improve SOD through relieving immune and inflammation dysfunction.

Additionally, we discovered three genes associated with the B cell receptor signaling pathway and complement system in the SYGC treatment of SOD, namely AURKB, KIF11 and PLG. Aurora kinase B (AURKB), which belongs to the mitotic protein kinase family, played a role in mitosis and the inflammatory pathway through of NF-κB transcription ^[39,40]. The plasminogen protein encoded by PLG can regulate skeletal muscle regeneration ^[41] and the resolution of inflammation through macrophage polarization and efferocytosis ^[42]. Moreover, the expression of these 3 genes returned to normal after SYGC treatment. These data imply that the three genes may play a critical role in the SYGC treatment of SOD.

Furthermore, glycycoumarin, licoflavonol, echinatin and homobutein exert good binding activity with the above three genes. It was reported that glycycoumarin can relax gastrointestinal smooth muscle tone ^[43] and inhibit tetanic contractions ^[38]. Echinatin and homobutein exerted favorable pharmacological effects on anti-inflammatory and anti-oxidant activity partly to NF-κB inhibition ^[44,45,46]. Collectively, these chemical components of SYGC provide the pharmacological basis for the immune and inflammation activities related to SOD.

5.   Conclusions

In conclusion, the present study displayed that multi-ingredient therapeutics of SYGC regulated SOD development by multi-targets and multi-pathways. Future studies should be conducted to explore the involvement of these targets in the treatment of SOD with SYGC. In addition, more experiments are still needed to confirm our findings. Nevertheless, our research still provided some reasonable major mediators for the anti-SOD effects of SYGC, including four active compounds (glycycoumarin, licoflavonol, echinatin and homobutein), three targets (AURKB, KIF11 and PLG) and several pathways (B cell receptor signaling pathway, complement system and muscle contraction).

Acknowledgments

This work was supported by grants from the National Natural Science Foundation of China (No. 81904017).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.