In this paper, a time periodic and two–group reaction–diffusion epidemic model with distributed delay is proposed and investigated. We firstly introduce the basic reproduction number for the model via the next generation operator method. We then establish the threshold dynamics of the model in terms of , that is, the disease is uniformly persistent if , while the disease goes to extinction if . Finally, we study the global dynamics for the model in a special case when all the coefficients are independent of spatio–temporal variables.
Citation: Lin Zhao, Zhi-Cheng Wang, Liang Zhang. Threshold dynamics of a time periodic and two–group epidemic model with distributed delay[J]. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1535-1563. doi: 10.3934/mbe.2017080
Related Papers:
[1]
Mario Grassi, Giuseppe Pontrelli, Luciano Teresi, Gabriele Grassi, Lorenzo Comel, Alessio Ferluga, Luigi Galasso .
Novel design of drug delivery in stented arteries:
A numerical comparative study. Mathematical Biosciences and Engineering, 2009, 6(3): 493-508.
doi: 10.3934/mbe.2009.6.493
[2]
Adam Peddle, William Lee, Tuoi Vo .
Modelling chemistry and biology after implantation of a drug-eluting stent. Part Ⅱ: Cell proliferation. Mathematical Biosciences and Engineering, 2018, 15(5): 1117-1135.
doi: 10.3934/mbe.2018050
[3]
Nilay Mondal, Koyel Chakravarty, D. C. Dalal .
A mathematical model of drug dynamics in an electroporated tissue. Mathematical Biosciences and Engineering, 2021, 18(6): 8641-8660.
doi: 10.3934/mbe.2021428
[4]
Xiaobing Feng, Tingao Jiang .
Mathematical and numerical analysis for PDE systems modeling intravascular drug release from arterial stents and transport in arterial tissue. Mathematical Biosciences and Engineering, 2024, 21(4): 5634-5657.
doi: 10.3934/mbe.2024248
[5]
Cameron Meaney, Gibin G Powathil, Ala Yaromina, Ludwig J Dubois, Philippe Lambin, Mohammad Kohandel .
Role of hypoxia-activated prodrugs in combination with radiation therapy: An in silico approach. Mathematical Biosciences and Engineering, 2019, 16(6): 6257-6273.
doi: 10.3934/mbe.2019312
[6]
Saied M. Abd El-Atty, Konstantinos A. Lizos, Osama Alfarraj, Faird Shawki .
Internet of Bio Nano Things-based FRET nanocommunications for eHealth. Mathematical Biosciences and Engineering, 2023, 20(5): 9246-9267.
doi: 10.3934/mbe.2023405
[7]
Hatim Machrafi .
Nanomedicine by extended non-equilibrium thermodynamics: cell membrane diffusion and scaffold medication release. Mathematical Biosciences and Engineering, 2019, 16(4): 1949-1965.
doi: 10.3934/mbe.2019095
[8]
Avner Friedman, Najat Ziyadi, Khalid Boushaba .
A model of drug resistance with infection by health care workers. Mathematical Biosciences and Engineering, 2010, 7(4): 779-792.
doi: 10.3934/mbe.2010.7.779
[9]
John Boscoh H. Njagarah, Farai Nyabadza .
Modelling the role of drug barons on the prevalence of drug epidemics. Mathematical Biosciences and Engineering, 2013, 10(3): 843-860.
doi: 10.3934/mbe.2013.10.843
[10]
R. A. Everett, Y. Zhao, K. B. Flores, Yang Kuang .
Data and implication based comparison of two chronic myeloid leukemia models. Mathematical Biosciences and Engineering, 2013, 10(5&6): 1501-1518.
doi: 10.3934/mbe.2013.10.1501
Abstract
In this paper, a time periodic and two–group reaction–diffusion epidemic model with distributed delay is proposed and investigated. We firstly introduce the basic reproduction number for the model via the next generation operator method. We then establish the threshold dynamics of the model in terms of , that is, the disease is uniformly persistent if , while the disease goes to extinction if . Finally, we study the global dynamics for the model in a special case when all the coefficients are independent of spatio–temporal variables.
1.
Introduction
The tumor microenvironment (TME) is a complex ecosystem comprising numerous cells, including immune cells, stromal cells and immunosuppressive cells [1]. Tumour endothelial cells (TECs) are a prominent component of TME and TECs, which have critical functions on tumor growth, progression and metastasis [2]. As tumor angiogenesis is essential for tumor growth and metastasis, inducing tumor-associated angiogenesis is a promising tactic in cancer progression [3]. Angiogenesis is now accepted as an important hallmark of cancer [4] and endothelial cells (which form tight adhesions to ensure vessel integrity) are essentials in inducing angiogenesis in TME [3]. In TME, TECs interact with tumor cells via juxtacrine and paracrine signaling during tumor intravasation and metastasis [5].
Colorectal cancer (CRC) is the fourth most common cancer and a leading cause of cancer mortality worldwide [6]. Recently stated that distinct stromal transcriptional and interactions are altered in colon cancer development [7] and stromal cells contribute to the CRC transcriptome [8]. It was also stated that stromal gene expression defines poor prognosis in colorectal cancer [9] and induces poor-prognosis subtypes in colorectal cancer [10]. Altogether, stromal cells have substantial influence and regulatory roles in colon cancer and endothelial cells (ECs) may cause colon carcinogenesis as a part of the stromal component.
So, analyzing the transcriptomes of colon TECs (CTECs) in TME reveals new targets and avenues and explores more uncontrolled factors in colon cancer research. We present a comprehensive analysis through bioinformatic tools and identified molecular and genetic alterations in TECs. We identified DEGs from microarray datasets of gene expression omnibus (GEO), hub genes from the interactions of deregulated genes and find out significant KEGG pathways. We also identified significant hub genes that are linked with poor survival of patients. Moreover, we found that hub genes are associated with immune cells infiltration and positively correlated with immune suppressive markers. This integrated analysis provides vital molecular insights into CTECs characterization, which may directly affect treatment recommendations for colon cancer patients. It provides opportunities for genome-guided clinical trials as well as drug development for the treating of colon cancer patients.
2.
Materials and methods
2.1. Datasets
We searched the NCBI gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) using the keywords "colon cancer endothelial cells, " "endothelial cells, " and "normal endothelial cells, " and identified one CTECs gene expression dataset GSE89287 (n = 24) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89287) [11]. This dataset included different cell types (endothelial cells, macrophages and epithelial cells) which were isolated from non-paired primary normal colon tissues and colorectal carcinomas and subsequently, RNA was isolated. In this study, we included only endothelial cells and excluded the other cell types from our study. Endothelial cells were isolated using the facs-sorting method. The dataset included 24 samples of endothelial cells, including 10 normal colon samples and 14 colon tumor samples. The platform of this data is GPL4133 (Agilent-014850 whole human genome microarray 4x44K G4112F) with a feature number version. Gene expression profiling interactive analysis (GEPIA) (http://gepia2.cancer-pku.cn/#index) dataset (n = 316) was used for the study of prognostic and expression levels of hub genes [12]. In addition, the correlation between the expression of hub genes and tumor-infiltrating immune cells was explored via gene modules in the TIMER dataset (n = 457) (https://cistrome.shinyapps.io/timer/) [13]. For identifying the gene-gene correlation, we used the TCGA-COAD dataset (n = 287) (https://gdc-portal.nci.nih.gov/) which was normalized into base log2.
2.2. Identification of differentially expressed genes (DEGs) between CTECs and CNECs
We used the web tool Network Analyst [14]to identify the significant DEGs between CTECs and CNECs. Dataset was normalized by quantile normalization and the R package "limma" was utilized to identify the DEGs between CTECs and CNECs. We utilized the Network Analyst tool to identify the average expression level of single gene having multiple probes in this gene expression study. We identified the DEGs with a threshold of |LogFC| > 0.585 and adjusted P value < 0.05.
2.3. Gene-set enrichment analysis
We performed gene-set enrichment analysis of the DEGs by GSEA [15]. We identified the KEGG [16] pathways that were significantly associated with the up-regulated and the down-regulated DEGs, respectively. We identified the significant pathways with a threshold of false discovery rate (FDR) < 0.05.
2.4. Construction of protein-protein interaction (PPI) network and modular analysis of PPI
To better know the relationship among these screened DEGs, the PPI network was established using the STRING database [17]. The hub genes in the PPI network were identified according to a degree using the node explorer module of Network Analyst software [14]. A hub gene was defined as a gene that was connected to no less than 10 other DEGs. A cytoscape plug-in molecular complex detection (MCODE) was employed to detect the modules from the PPI network [18]. We identified the significant modules based on the MCODE score and node number. The threshold of the MCODE was Node Score Cut-off: 0.2, Haircut: true, K-Core: 2 and maximum depth from Seed: 100.
2.5. Survival analysis of hub genes
We compared the overall survival (OS) and the disease-free survival (DFS) of colon cancer patients classified based on gene expression levels (expression levels > median versus expression levels < median). Kaplan-Meier survival curves were used to show the survival differences and the log-rank test was utilized to evaluate the significance of survival differences. The prognostic roles of screened hub genes in COAD were analyzed using expression profiling interactive analysis (GEPIA) databases [12]. Cox P < 0.05 was considered as significant between two groups of patients.
2.6. Validation of prognosis-associated hub gene expression levels in TCGA COAD dataset
The expression levels of top hub genes were further validated using the gene expression profiling interactive analysis (GEPIA), a newly developed interactive web server for analyzing the RNA sequencing data [12]. We used TCGA COAD tumor data with matched normal. Hub genes with |log2FC| > 0.585 and P < 0.05 were considered statistically significant between the two groups.
2.7. Evaluation of immune scores and stromal scores in colon cancer
We utilized the "ESTIMATE" R package to quantify the immune scores (predict the immune cell quantity) and stromal scores (predict the stromal cells quantity) for each of the tumor samples [19]. Then we calculated Spearman's correlations between the expression levels of prognostic hub genes and immune and stromal scores. The threshold value of correlation is R > 0.30 and P-value is less than 0.001 (Spearman's correlation test).
2.8. Analysis of immune cell infiltration and correlation of top hub genes with immune-inhibitory markers
We analyzed the significant correlations of hub genes with the abundance of immune infiltrates, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells, via gene modules in TIMER [13]. In addition, Spearman's correlations between the expression level of hub genes and the immune inhibitory marker genes of tumor-infiltrating immune cells were explored via correlation modules in TIMER [13]. Moreover, we identified pearson's correlations between the expression level of hub genes and the marker genes of inhibitory immune cells in TCGA-COAD datasets using R software. The R package "ggplot2" was employed for drawing a correlation graph [20]. The gene markers of tumor-infiltrating immune cells included markers of monocytes, TAMs, M2 macrophages, T-helper 1 (Th1) cells, Tregs and exhausted T cells. These gene markers are illustrated in prior studies [12,21,22,23,24].
2.9. Statistical analysis
We used R programming software for the statistical analysis of this study. For the differential expression analysis of gene expression data, only genes with |logFC| > 0.585 and adjusted P-value < 0.05 were considered statistically significant between the two groups. In the survival analysis, Cox P < 0.05 was considered as statistically significant. Spearman's correlation test between the ssGSEA scores and the expression level of specific prognostic genes was performed because these data were not normally distributed (P < 0.05) [25]. For analyzing the correlations between the expression levels of hub genes with the expression levels of other marker genes, we employed Pearson's correlation test because these data were normally distributed [25].
3.
Results
3.1. Identification of Differentially Expressed Genes
Based on the log2FC and adjusted P-value, we identified 362 differentially expressed genes (DEGs) in CTECs relative to the CNECs, which included 117 up-regulated genes (Supplementary Table S1) and 245 down-regulated (Supplementary Table S2) in the CTECs samples. Table 1 demonstrates the top ten up-regulated genes (SPARC, IGFBP7, COL4A2, SELE, CXCL1, VWF, CCDC3, FSTL1, VWA1 and IFITM3) with significant statistical description. The most up-regulated gene SPARC is up-expressed in the stromal portion of CRC tissues [26]. The extracellular matrix-associated gene, COL4A2, is linked with cancer stemness [27]. Another up-regulated gene, CXCL1, is associated with increasing the metastatic ability of colon cancer by enhancing cell migration, matrix metalloproteinase-7 expression and EMT [28]. VWF, another up-regulated gene, is a plasma marker for the early detection of adenoma and colon cancer [29]. Table 2 illustrated the top ten down-regulated genes (IGKV1-5, IGLV1-44, IGKC, IGHV3-30, KRT81, IGHG1, IGHM, POU2AF1, IGKV2-24 and IGLJ3) with all significant statistical descriptions. Interestingly, most of the top down-regulated genes are associated with immunological activities. Some other down-regulated genes of CTECs have been demonstrated to be under-expressed in CRC. For example, the expression of POU2AF1 is down-regulated in the rat colon and in turn, reveals proximal-distal differences in the process of histone modifications and proto-oncogene expression [30]. Altogether, abnormally expressed numerous genes in CTECs compared to NTECs identified by the bioinformatic analysis have been associated with CRC pathology and carcinogenesis.
3.2. KEGG pathway enrichment analysis of significantly identified DEGs
The functional enrichment analysis identifying the enriched up-regulated and down-regulated biological pathways that are associated with significant DEGs (Figure 1). GSEA pathway analysis revealed nine up-regulated pathways (focal adhesion, ECM-receptor interaction, pathways in cancer, small cell lung cancer, adherens junction, arrhythmogenic right ventricular cardiomyopathy (ARVC), pathogenic escherichia coli infection, cytokine-cytokine receptor interaction and proximal tubule bicarbonate reclamation) (Figure 1A) in CTECs. In addition, we also identified ten down-regulated pathways (protein export, cell adhesion molecules (CAMs), arrhythmogenic right ventricular cardiomyopathy (ARVC), N-Glycan biosynthesis, the intestinal immune network for IgA production, primary immunodeficiency, hypertrophic cardiomyopathy (HCM), ECM-receptor interaction, dilated cardiomyopathy and vasopressin-regulated water reabsorption) in CTECs (Figure 1B). Recently, it was stated that focal adhesion, ECM-receptor interaction, pathways in cancer, small cell lung cancer, adherens junction and cytokine-cytokine receptor interaction pathways are dysregulated in colon tumors stroma [7].
Figure 1.
KEGG pathway enrichment analysis of significant DEGs. A. Significantly up-regulated pathways in colon tumor endothelial cells. B. Significantly down-regulated pathways in colon tumor endothelial cells. FDR: false discovery rate.
3.3. Construction of PPI by using significant DEGs and identification of hub genes and functional clusters from the PPI
Using all 362 DEGs, PPI networks were constructed using STRING software and hub genes were identified using the node explorer of network analyst software. Many significant DEGs are involved in the networks (Supplementary Table S3). The identification of hub genes explored that FN1, GAPDH, COL1A1, COL1A2, CTNNB1, LAMB1, LAMC1, SPARC, IGFBP3 andCOL4A1 are top up-regulated hub genes and HSP90B1, PDIA6, ITGA4, DDOST, SDC1, CKAP4, CD38, MANF, ITGB7 andCCR2 are top down-regulated hub genes in CTECs respectively (Figure 2 and Table 3). Interestingly, we found that top-up-regulated (Table 1) SPARC also acts as a hub gene in the PPI network. FN1 suppressed apoptosis and promoted viability, invasion and migration of tumor cells in CRC [31]. It was showed that COL1A1 and COL1A2 are associated with colon cancer [32,33]. The high expression levels of CTNNB1 positively correlated with metastasis of colon cancer [34]. Deregulation of SDC1 contributes to cancer progression by promoting cell proliferation, metastasis, invasion and angiogenesis [35]. Altogether, it may be stated that CTECs-derived these hub genes may contribute to colon cancer pathogenesis.
Figure 2.
The functional enrichment analysis of MCODE derived clusters. A. PPI interaction of Cluster 1. B. Enrichment of the KEGG pathways in Cluster 1. C. PPI interaction of Cluster 3. D. Enrichment of the 11 KEGG pathways in Cluster 3. E. Interaction of 23 nodes in Cluster 4. F. The enrichment of 11 KEGG pathways in Cluster 4.
Moreover, MCODE identified 15 clusters from the original PPI networks. The description of MCODE derived clusters is illustrated in Table 4. The top significant cluster 1 contained 15 nodes and 105 edges (Figure 2 and Table 4). We identified the functional enrichment of KEGG pathways for all clusters by using the GSEA. Interestingly, we found that eight of the clusters out of 15 are associated with the enrichment of KEGG pathways. Gene set of Cluster 1 is related to the enrichment of 4 pathways: ECM-receptor interaction, small cell lung cancer, pathways in cancer and focal adhesion (Figure 2A, B). Cluster 3 (Figure 2C, D), Cluster 4 (Figure 2E, F) and Cluster 9 (Cell adhesion molecules and natural killer cell-mediated cytotoxicity) are associated with the enrichment of immune, stromal and cancer-associated pathways. Some of the enriched pathways in Cluster 3 included ECM-receptor interaction, focal adhesion, small cell lung cancer, cell adhesion molecules, the intestinal immune network for IgA production, pathways in cancer and hematopoietic cell lineage. In addition, some of the pathways enriched in Cluster 4 are the chemokine signaling pathway, adherens junction, cytokine-cytokine receptor interaction, focal adhesion, NOD-like receptor signaling pathway, pathways in cancer, leukocyte trans-endothelial migration and tight junction. Protein export is enriched in Cluster 2 and the B cell receptor signaling pathway is enriched in Cluster 5. Cluster 13 (Fatty acid metabolism and PPAR signaling pathway) and cluster 15 (Glycosphingolipid biosynthesis-Lacto and neolacto series and N-Glycan biosynthesis) are mainly associated with the enrichment of metabolic pathways. Uddin et al. also identified some of these pathways, including pathways in cancer, small cell lung cancer, ECM-receptor interaction, focal adhesion, natural killer cell-mediated cytotoxicity and cell adhesion molecules, are enriched in colon tumor stroma [7,9]. Altogether, it can be stated that this hub PPI from CTECs may contribute to colon carcinogenesis.
Table 4.
MCODE identified 15 clusters from the PPI networks.
3.4. Hub genes are negatively associated with prognosis
We considered the top ten up-regulated and top ten down-regulated genes from the original PPI network for survival analysis. We selected the TCGA-COAD database for screening OS and DFS. We got COL1A1, COL1A2, IGFBP3, SPARC andDDOST hub genes are significantly involved with patient survival time (Figure 3). We found that the higher expression level of up-regulated COL1A1, COL1A2, IGFBP3 and SPARC are worse for patient's DFS and COL1A2 is also associated with poor OS. In addition, the low expression level of down-regulated DDOST is more inferior for patient's OS and DFS. It showed that higher expressions of COL1A1 were related to poor survival in colon cancer patients [33]. COL1A2 gene is also associated with the prognosis of IIA stage colon cancer [32]. Stromal SPARC had a pro-metastatic impact in vitro and was a characteristic of aggressive tumors with a poor prognosis in CRC patients [26]. Collectively, hub genes of CTECs are prognostic markers in COAD and may be contributed to colon carcinogenesis.
Figure 3.
Survival analysis of individual hub genes. Up-regulated COL1A1, COL1A2, IGFBP3 and SPARC and down-regulated DDOST have been associated with the poor prognosis in COAD.
3.5. Comparisons of the expression levels of prognostic hub genes between colon cancer and normal tissue
Interestingly, we found that all four up-regulated hub genes (COL1A1, COL1A2, IGFBP3 and SPARC) were consistently up-regulated in TCGA COAD samples versus normal samples (P < 0.05) (Figure 4). It indicates that CTECs may have a contribution to the COAD-derived transcriptomes in colon cancer. However, one of the down-regulated hub genes, DDOST, showed no significant expression differences between TCGA-COAD samples versus normal samples (P < 0.05), suggesting that this hub gene is expressed explicitly in CTECs.
Figure 4.
Validation of the expression levels of prognostic hub genes between colon cancer and normal tissue (P < 0.05, |Log FC| > 0.585).
3.6. Prognostic hub genes are correlated with immune infiltrations and associated with immune-inhibitory markers in colon TME
Tumor-infiltrating lymphocytes are independent prognostic factors for better survival and sentinel lymph node status in cancers [36]. Genes highly expressed in the microenvironment are expected to negatively associate with tumor purity, while the opposite is expected for genes highly expressed in the tumor cells [13]. So, we assessed whether the expression of hub genes was correlated with immune infiltration levels in COAD. We investigated the correlations of prognostic hub genes (COL1A1, COL1A2, IGFBP3, SPARC and DDOST) with immune scores and stromal scores in COAD. Interestingly, we found that the expression levels of COL1A1, COL1A2, IGFBP3 and SPARC positively correlated with the immune scores (Figure 5A) and stromal scores (Figure 5B) (Spearman's correlation test, R > 0.30, P < 2.2e−16). This result indicated that the prognostic hub genes are associated with the regulation of the tumor microenvironment.
Figure 5.
Up-regulated prognostic hub genes are associated with the regulation of the tumor microenvironment in COAD. The expression level of COL1A1, COL1A2, IGFBP3 and SPARC positively correlated with the immune scores (A) and stromal scores (B).
In addition, we identified the correlations between the expression of hub genes and the immune infiltration levels in COAD by using the TIMER tool. The results discovered that up-regulated COL1A1, COL1A2, IGFBP3 and SPARC expression have significant (P ≤ 0.05) correlations with tumor purity in COAD but down-regulated DDOST not correlated with tumor purity (Figure 6). In addition, COL1A1 expression has significant correlations with infiltrating levels of CD4+ T cells, macrophage, neutrophil and dendritic cells. We also found that a significant positive correlation between COL1A2, IGFBP3 and SPARC expression and infiltration of all five immune cells (CD4+ T cells, CD8+ T cells, Neutrophils, Macrophages and Dendritic cells) in COAD. Besides, we found significant negative correlations between DDOST expression and infiltration of all five immune cells (CD4+ T cells, CD8+ T cells, Neutrophils, Macrophages and Dendritic cells) (Figure 6). Recently it was stated that macrophages and neutrophils are associated with immunosuppressive inflammation to modulate anti-tumor immunity [37]. DCs can promote tumor metastasis by increasing Treg cells and reducing CD8+ T cell cytotoxicity [38]. These results suggest that the expression of COL1A1, COL1A2, IGFBP3, SPARC and DDOST has potential roles in tumor purity and immune cell infiltrations in colon cancer.
Figure 6.
Hub genes are correlated with immune infiltration level and tumor purity. Expression of COL1A1, COL1A2, SPARC, IGFBP3 and DDOST are significantly negatively correlated to tumor purity and have significant positive correlations with infiltrating levels of CD4+ T cells, macrophages, neutrophils and dendritic cells in COAD. The correlation analyses were performed using TIMER [13]. RSEM: RNA-Seq by Expectation Maximization [39].
To elucidate the relationship between up-regulated hub genes COL1A1, COL1A2, IGFBP3 and SPARC and the diverse immune infiltrating cells, we find out the correlations between these genes and immune marker sets of various immune cells including monocytes, TAMs, M2 macrophages, Th1, Tregs and T cell exhaustion (Table 5). Surprisingly, we got that the hub genes COL1A1, COL1A2, IGFBP3 andSPARC are positively correlated (Pearson correlation test, P ≤ 0.001) with the immune markers of monocytes, TAMs, M2 macrophages, Th1, Tregs and T cell exhaustion. These findings suggest that the high expression level of hub genes plays an essential role in the infiltration of monocytes, TAMs, M2 macrophages, Th1, Tregs and T cell exhaustion in COAD. It was recently reported that LAYN acts as a prognostic biomarker for determining prognosis and immune infiltration in gastric and colon cancers [21]. Potent immunosuppressive T-regulatory cells (Tregs) are found in a vast array of tumor types and tumor-infiltrating Tregs are often associated with poor clinical outcomes. Tregs also promote cancer progression through their ability to limit anti-tumor immunity and promote angiogenesis [40]. Another marker, FOXP3, plays an important role in Treg cells which leads to the suppression of cytotoxic T cells attacking tumor cells [40]. TIM-3 is a crucial surface protein on exhausted T cells, a critical gene that regulating T cell exhaustion [41].
Table 5.
Correlation analysis between survival-associated hub genes and markers of immune-inhibitory cells in TCGA-COAD data.
Immune cell
Gene markers
SPARC
COL1A1
COL1A2
IGFBP3
R
P
R
P
R
P
R
P
Monocyte
CD86
0.69
1.49E−42
0.65
2.71E−35
0.67
3.02E−39
0.55
2.39E−24
CSF1R
0.71
1.22E−44
0.68
5.80E−41
0.72
1.60E−46
0.59
9.29E−28
TAM
CCL2
0.76
3.55E−56
0.69
1.00E−41
0.72
2.51E−47
0.58
2.49E−27
CD68
0.49
8.62E−19
0.48
5.34E−18
0.49
1.16E−18
0.46
1.25E−16
IL10
0.57
1.86E−26
0.52
4.70E−21
0.54
2.65E−23
0.42
1.87E−13
M2 macrophage
CD163
0.71
4.72E−46
0.71
4.11E−46
0.73
7.07E−49
0.52
1.88E−21
VSIG4
0.73
8.85E−50
0.70
2.16E−43
0.71
6.58E−46
0.56
7.05E−25
MS4A4A
0.70
2.77E−44
0.65
1.77E−35
0.67
7.62E−39
0.54
6.22E−23
Th1
T-bet (TBX21)
0.36
3.16E−10
0.38
1.66E−11
0.38
3.46E−11
0.37
8.76E−11
IFN-γ (IFNG)
0.18
0.001
0.21
0.00040
0.20
0.000834
0.16
0.00589
TNF-α (TNF)
0.40
3.37E−12
0.38
1.65E−11
0.40
3.16E−12
0.21
0.00030
Treg
FOXP3
0.56
6.11E−25
0.55
3.52E−24
0.58
1.49E−27
0.50
1.01E−19
CCR8
0.56
1.13E−24
0.54
4.14E−23
0.58
2.02E−27
0.48
1.23E−17
TGFβ (TGFB1)
0.74
1.08E−50
0.73
1.99E−49
0.75
1.61E−52
0.62
6.64E−32
T cell exhaustion
PD-1 (PDCD1)
0.32
3.27E−08
0.32
2.35E−08
0.33
9.61E−09
0.35
8.10E−10
CTLA4
0.42
1.40E−13
0.42
1.10E−13
0.44
7.38E−15
0.40
3.58E−12
LAG3
0.25
1.56E−05
0.29
6.10E−07
0.28
1.94E−06
0.27
3.07E−06
TIM-3(HAVCR2)
0.70
1.48E−43
0.67
2.34E−38
0.69
3.05E−41
0.54
1.58E−23
TIGIT
0.39
4.84E−12
0.41
7.99E−13
0.42
5.46E−14
0.40
2.49E−12
CXCL13
0.35
1.82E−09
0.36
5.46E−10
0.35
8.72E−10
0.31
1.34E−07
LAYN
0.89
2.42E−97
0.79
1.39E−63
0.83
4.05E−75
0.67
2.35E−38
*Note: R is Pearson correlation and P is p-value in Pearson's correlation test.
Since the elevated expression of prognostic hub genes were positively associated with the immune inhibitory markers PD-1 and TGFB1 (Table 5), we expect the expression levels of PD-L1 (CD274), PD-L2 (PDCD1LG2) and TGFBR1 could be positively associated with the prognosis-associated hub genes. Interestingly, we found that PDL1 (CD274), PDL2 (PDCD1LG2) and TGFBR1 are moderate to highly correlated with all four-prognosis related highly expressed hub genes: COL1A1, COL1A2, SPARC and IGFBP3 (Figure 7). T cell immunity recovered by PD-1/PD-L1 immune checkpoint blockade has been demonstrated to be a promising cancer therapeutic strategy [42]. The expression of PD-L2 was observed in tumors, stroma and endothelial cells and PD-L2 expression is relevant to anti-PD-1 therapy in cancer [43]. TGFBR1 is one of the receptors for TGF-β ligands. The TGF-β signaling pathway is crucially associated with stimulation, induction and maintenance of EMT in cancers [44]. TGFBR1 is up-regulated in multiple malignancies and dysregulation of TGFBR1 leads to tumorigenesis by controlling cellular signaling, which in turn is associated with colorectal cancer risk [45]. Altogether, it suggests that the prognostic associated hub genes may regulate the immune-suppressive activities of TME through the interactions with the immune marker in colon cancer.
Figure 7.
Prognostic hub genes are significantly positively correlated with immune suppressive PD-L1, PD-L2 and TGFBR1. Up-regulated four hub genes (SPARC, COL1A2, COL1A2 and IGFBP3) are significantly correlated with three prominent immune-inhibitory markers (PD-L1, PD-L2 and TGFBR1). The correlation analyses were performed using R software. The Pearson's correlation test p-values (P) and correlation coefficients (R) are shown in the figure.
TECs are the substantial component of tumor stroma in TME and these cells are associated with tumor malignancies, progression and metastasis [5]. This present study conducted a differential expression analysis using gene expression profiling data from the GEO database. We identified 362 DEGs in CTECs, including 117 up-regulated genes (SPARC, COL4A2, CXCL1 andVWF) and 245 down-regulated genes. Interestingly, we found that many immunological genes are down-regulated in CTECs (Table 2). Subsequent pathway enrichment analysis revealed that up-regulated DEGs were significantly enriched with cancer, cellular development and immune regulation. Down-regulated DEGs were enriched considerably with mostly immunodeficiency and metabolism-related pathways (Figure 1). These results revealed the abnormal cellular growth, immune regulatory, cancerous and metabolism-associated pathways in CTECs.
Next, we employed DEGs to construct a PPI network and extracted significant clusters from the original PPI network (Figure 2). Finally, we identified hub genes (degree > 10) that are dysregulated in CTECs. We identified CTECs-derived five hub genes (COL1A1, COL1A2, SPARC, IGFBP3 andDDOST) whose expression was significantly associated with the poor prognosis in TCGA-COAD data (Figure 3). These hub genes are mainly involved in protein digestion and absorption (COL1A1 and COL1A2), cellular signaling (SPARC and IGFBP3) and enzymatic action of metabolism (DDOST), suggesting that the deregulation of these cellular functions in CTECs may contribute to the altered prognosis in colon cancer. In addition, the expression analysis of these prognostic hub genes in colon cancer was further evaluated using the GEPIA database [12]. We found that four hub genes (COL1A1, COL1A2, SPARC and IGFBP3) are also significantly up-regulated in TCGA-COAD samples. Still, another down-regulated hub gene (DDOST) is not considerably down-regulated in TCGA-COAD samples (Figure 4). It indicates that CTECs-associated gene signature substantially contributed to colon carcinogenesis.
Furthermore, our study showed that the immune infiltration levels and diverse immune marker sets are correlated with the expression levels of COL1A1, COL1A2, SPARC and IGFBP3 hub genes. Our analysis demonstrated that a significant positive correlation between the COL1A1, COL1A2, SPARC and IGFBP3 expression levels and infiltration level of CD4+ T cells, macrophages, neutrophils and DCs in COAD (Figure 6). These results are suggesting that the potential immunological roles of hub genes in colon TME. In addition, our study indicated that COL1A1, COL1A2, SPARC and IGFBP3 could activate Tregs and induce T cell exhaustion. The increase in up-regulated four hub genes expression positively correlates with Treg and T cell exhaustion markers (Table 5). Other markers of monocytes, TAM, M2 macrophage, Th1, are also positively correlated with the expression of these hub genes (Table 5). Some other crucial immune inhibitory markers, PD-L1, PD-L2 and TGFBR1, are significantly correlated with the expression level of COL1A1, COL1A2, SPARC and IGFBP3 (Figure 7). Altogether, the prognostic hub genes are associated with immunosuppressive colon TME.
Since our analysis identified numerous CTECs-derived transcriptomes that are critically associated with colon cancer, we speculated that the endothelial cells might contribute to the other vascular pathologies, including the central nervous system (CNS) and other diseases. It was stated that cerebral endothelial cells play an active role in the pathogenesis of CNS inflammatory diseases [46]. Endothelial cells are associated with the leukocyte migration in CNS and regulating the leukocyte-endothelial cell interactions and the crosstalk between endothelial cells and glial cells or platelets in CNS [46]. The major vascular anomalies of the nervous system included cerebral cavernous malformations (CCMs) and arteriovenous malformations (AVMs) [47,48]. These CNS-associated vascular diseases are characterized by genetic alterations, such as gene polymorphisms and mutations [47,48,49,50,51]. For example, germline mutation enrichment in pathways controlling the endothelial cell homeostasis in patients with brain arteriovenous malformation [48]. Altogether, our findings may provide insights to find the roles of endothelial cells in diseases, including cancer, CNS-associated disease and other vascular pathologies.
This study identified numerous deregulated transcriptomes in the CTECs that could be used as biomarkers for the diagnosis and prognosis of colon cancer and may provide therapeutic targets for colon cancer. However, to translate these findings into clinical application, further experimental and clinical validation would be necessary.
5.
Conclusions
In summary, we identified potential CTECs-derived significantly deregulated transcriptomes that are involving with the pathogenesis, prognosis and immune inhibition within the tumor microenvironment of colon cancer. This study reveals a potential regulatory mechanism of CTECs in the tumor microenvironment of colon cancer and may contribute to revealing CTECs-colon cancer cellular cross-talk.
Conflicts of interest
The authors declare no conflict of interest.
Acknowledgments
We thank all laboratory staff and technicians for their participation in this study. The present study was supported by the project from Xinjiang Health Poverty Alleviation Work Appropriate Technology Promotion (Grant no. SYTG-Y202035).
Author contributions
Jie Wang and Md. Nazim Uddin contributed equally to this work. Jie Wang and Md. Nazim Uddin conceived the research, designed analysis strategies. Rehana Akter wrote the manuscript and Yun Wu conceived the study and wrote the manuscript. All the authors read and approved the final manuscript.
References
[1]
[ S. Altizer,A. Dobson,P. Hosseini,P. Hudson,M. Pascual,P. Rohani, Seasonality and the dynamics of infectious disease, Ecol. Lett., 9 (2006): 467-484.
[2]
[ R. M. Anderson, Discussion: the Kermack-McKendrick epidemic threshold theorem, Bull. Math. Biol., 53 (1991): 3-32.
[3]
[ R. M. Anderson and R. May,
Infectious Diseases of Humanns: Dynamics and Control, Oxford University Press, Oxford, 1991.
[4]
[ N. Bacaër,D. Ait,H. El, Genealogy with seasonality, the basic reproduction number, and the influenza pandemic, J. Math. Biol., 62 (2011): 741-762.
[5]
[ N. Bacaër,S. Guernaoui, The epidemic threshold of vector–borne disease with seasonality, J. Math. Biol., 53 (2006): 421-436.
[6]
[ E. Beretta,T. Hara,W. Ma,Y. Takeuchi, Global asymptotic stability of an SIR epidemic model with distributed time delay, Nonlinear Anal., 47 (2001): 4107-4115.
[7]
[ B. Bonzi,A. A. Fall,A. Iggidr,G. Sallet, Stability of differential susceptibility and infectivity epidemic models, J. Math. Biol., 62 (2011): 39-64.
[8]
[ F. Brauer, Compartmental models in epidemiology, Mathematical Epidemiology, Springer, 56 (2008): 19-79.
[9]
[ L. Burlando, Monotonicity of spectral radius for positive operators on ordered Banach space, Arch. Math., 56 (1991): 49-57.
[10]
[ L. Cai,M. Martcheva,X.-Z. Li, Competitive exclusion in a vector-host epidemic model with distributed delay, J. Biol. Dyn., 7 (2013): 47-67.
[11]
[ D. Dancer and P. Koch Medina, Abstract ecolution equations, Periodic problem and applications, Longman, Harlow, UK, 1992.
[12]
[ O. Diekmann,J. Heesterbeek,J. Metz, On the definition and the computation of the basic reproduction ratio in models for infectious diseases in heterogeneous populations, J. Math. Biol., 28 (1990): 365-382.
[13]
[ W. E. Fitzgibbon,M. Langlais,M. E. Parrott,G. F. Webb, A diffusive system with age dependency modeling FIV, Nonlinear Anal., 25 (1995): 975-989.
[14]
[ W. E. Fitzgibbon,C. B. Martin,J. J. Morgan, A diffusive epidemic model with criss–cross dynamics, J. Math. Anal. Appl., 184 (1994): 399-414.
[15]
[ W. E. Fitzgibbon,M. E. Parrott,G. F. Webb, Diffusion epidemic models with incubation and crisscross dynamics, Math. Biosci., 128 (1995): 131-155.
[16]
[ D. Gao and S. Ruan, Malaria models with spatial effects, John Wiley & Sons. (in press)
[17]
[ I. Gudelj,K. A. J. White,N. F. Britton, The effects of spatial movement and group interactions on disease dynamics of social animals, Bull. Math. Biol., 66 (2004): 91-108.
[18]
[ Z. Guo,F.-B. Wang,X. Zou, Threshold dynamics of an infective disease model with a fixed latent period and non–local infections, J. Math. Biol., 65 (2012): 1387-1410.
[19]
[ P. Hess, Periodic–Parabolic Boundary Value Problems and Positivity, Longman Scientific and
Technical, Harlow, UK, 1991.
[20]
[ H. Hethcote, The mathematics of infectious diseases, SIAM Rev., 42 (2000): 599-653.
[21]
[ W. Huang,K. Cooke,C. Castillo-Chavez, Stability and bifurcation for a multiple–group model for the dynamics of HIV/AIDS transmission, SIAM J. Appl. Math., 52 (1992): 835-854.
[22]
[ G. Huang,A. Liu, A note on global stability for a heroin epidemic model with distributed delay, Appl. Math. Lett., 26 (2013): 687-691.
[23]
[ J. M. Hyman,J. Li, Differential susceptibility epidemic models, J. Math. Biol., 50 (2005): 626-644.
[24]
[ H. Inaba, On a new perspective of the basic reproduction number in heterogeneous environments, J. Math. Biol., 65 (2012): 309-348.
[25]
[ Y. Jin,X.-Q. Zhao, Spatial dynamics of a nonlocal periodic reaction–diffusion model with stage structure, SIAM J. Math. Anal., 40 (2009): 2496-2516.
[26]
[ T. Kato,
Peturbation Theory for Linear Operators, Springer-Verlag, Berlin, Heidelerg, 1976.
[27]
[ J. Li,X. Zou, Generalization of the Kermack–McKendrick SIR model to a patchy environment for a disease with latency, Math. Model. Nat. Phenom., 4 (2009): 92-118.
[28]
[ J. Li,X. Zou, Dynamics of an epidemic model with non–local infections for diseases with latency over a patchy environment, J. Math. Biol., 60 (2010): 645-686.
[29]
[ M. Li,Z. Shuai,C. Wang, Global stability of multi–group epidemic models with distributed delays, J. Math. Anal. Appl., 361 (2010): 38-47.
[30]
[ X. Liang,X.-Q. Zhao, Asymptotic speeds of spread and traveling waves for monotone semiflows with applications, Comm. Pure Appl. Math., 60 (2007): 1-40.
[31]
[ Y. Lou,X.-Q. Zhao, Threshold dynamics in a time–delayed periodic SIS epidemic model, Discrete Contin. Dyn. Syst. Ser. B, 12 (2009): 169-186.
[32]
[ Y. Lou,X.-Q. Zhao, A reaction–diffusion malaria model with incubation period in the vector population, J. Math. Biol., 62 (2011): 543-568.
[33]
[ Y. Lou,X.-Q. Zhao, A theoretical approach to understanding population dynamics with deasonal developmental durations, J Nonlinear Sci., 27 (2017): 573-603.
[34]
[ P. Magal,C. McCluskey, Two–group infection age model including an application to nosocomial infection, SIAM J. Appl. Math., 73 (2013): 1058-1095.
[35]
[ P. Magal,X.-Q. Zhao, Global attractors and steady states for uniformly persistent dynamical systems, SIAM J. Math. Anal., 37 (2005): 251-275.
[36]
[ M. Martcheva,
An Introduction to Mathematical Epidemiology, Texts in Applied Mathematics, Springer, New York, 2015.
[37]
[ R. Martain,H. L. Smith, Abstract functional differential equations and reaction–diffusion system, Trans. Amer. Math. Soc., 321 (1990): 1-44.
[38]
[ C. McCluskey, Complete global stability for an SIR epidemic model with delay-distributed or discrete, Nonlinear Anal. Real World Appl., 11 (2010): 55-59.
[39]
[ C. McCluskey,Y. Yang, Global stability of a diffusive virus dynamics model with general incidence function and time delay, Nonlinear Anal. Real World Appl., 25 (2015): 64-78.
[40]
[ J. D. Murray,
Mathematical Biology, Springer-Verlag, Berlin, 1989.
[41]
[ R. Peng,X.-Q. Zhao, A reaction–diffusion SIS epidemic model in a time–periodic environment, Nonlinearity, 25 (2012): 1451-1471.
[42]
[ B. Perthame,
Parabolic Equations in Biology, Springer, Cham, 2015.
[43]
[ L. Rass and J. Radcliffe,
Spatial Deterministic Epidemics, Mathematical Surveys and Monographs, 102. American Mathematical Society, Providence, RI, 2003.
[44]
[ R. Ross, An application of the theory of probabilities to the study of a priori pathometry: Ⅰ, Proc. R. Soc. Lond., 92 (1916): 204-230.
[45]
[ S. Ruan, Spatial−temporal dynamics in nonlocal epidemiological models, Mathematics for Life Science and Medicine, Springer−Verlag, Berlin, (2007), 99–122.
[46]
[ S. Ruan and J. Wu, Modeling Spatial Spread of Communicable Diseases Involving Animal Hosts, Chapman & Hall/CRC, Boca Raton, FL, (2009), 293–316.
[47]
[ H. L. Smith,
Monotone Dynamical System: An Introduction to the Theorey of Competitive and Cooperative Systems, Math. Surveys and Monogr. vol 41, American Mathematical Society, Providence, 1995.
[48]
[ R. Sun, Global stability of the endemic equilibrium of multigroup SIR models with nonlinear incidence, Comput. Math. Appl., 60 (2010): 2286-2291.
[49]
[ Y. Takeuchi,W. Ma,E. Beretta, Global asymptotic properties of a delay SIR epidemic model with finite incubation times, Nonlinear Anal., 42 (2000): 931-947.
[50]
[ H. R. Thieme, Mathematics in population biology, Princeton University Press, Princeton, NJ, 2003.
[51]
[ H. R. Thieme, Spectral bound and reproduction number for infinite–dimensional population structure and time heterogeneity, SIAM J. Appl. Math., 70 (2009): 188-211.
[52]
[ H. R. Thieme, Global stability of the endemic equilibrium in infinite dimension: Lyapunov functions and positive operators, J. Differential Equations, 250 (2011): 3772-3801.
[53]
[ P. van den Driessche,X. Zou, Modeling relapse in infectious diseases, Math. Biosci., 207 (2007): 89-103.
[54]
[ B.-G. Wang,W.-T. Li,Z.-C. Wang, A reaction–diffusion SIS epidemic model in an almost periodic environment, Z. Angew. Math. Phys., 66 (2015): 3085-3108.
[55]
[ B.-G. Wang,X.-Q. Zhao, Basic reproduction ratios for almost periodic compartmental epidemic models, J. Dynam. Differential Equations, 25 (2013): 535-562.
[56]
[ L. Wang,Z. Liu,X. Zhang, Global dynamics of an SVEIR epidemic model with distributed delay and nonlinear incidence, Appl. Math. Comput., 284 (2016): 47-65.
[57]
[ W. Wang,X.-Q. Zhao, Threshold dynamics for compartmental epidemic models in periodic environments, J. Dynam. Differential Equations, 20 (2008): 699-717.
[58]
[ W. Wang,X.-Q. Zhao, A nonlocal and time-delayed reaction-diffusion model of dengue transmission, SIAM J. Appl. Math., 71 (2011): 147-168.
[59]
[ W. Wang,X.-Q. Zhao, Basic reproduction numbers for reaction–diffusion epidemic models, SIAM J. Appl. Dyn. Syst., 11 (2012): 1652-1673.
[60]
[ W. Wang,X.-Q. Zhao, Spatial invasion threshold of Lyme disease, SIAM J. Appl. Math., 75 (2015): 1142-1170.
[ D. Xu,X.-Q. Zhao, Dynamics in a periodic competitive model with stage structure, J. Math. Anal. Appl., 311 (2005): 417-438.
[63]
[ Z. Xu,X.-Q. Zhao, A vector–bias malaria model with incubation period and diffusion, Discrete Contin. Dyn. Syst. Ser. B, 17 (2012): 2615-2634.
[64]
[ L. Zhang,J.-W. Sun, Global stability of a nonlocal epidemic model with delay, Taiwanese J. Math., 20 (2016): 577-587.
[65]
[ L. Zhang and Z. -C. Wang, A time-periodic reaction-diffusion epidemic model with infection period, Z. Angew. Math. Phys. , 67 (2016), Art. 117, 14 pp.
[66]
[ L. Zhang,Z.-C. Wang,Y. Zhang, Dynamics of a reaction-diffusion waterborne pathogen model with direct and indirect transmission, Comput. Math. Appl., 72 (2016): 202-215.
[67]
[ L. Zhang,Z.-C. Wang,X.-Q. Zhao, Threshold dynamics of a time periodic reaction–diffusion epidemic model with latent period, J. Differential Equations, 258 (2015): 3011-3036.
[68]
[ Y. Zhang,X.-Q. Zhao, A reaction–diffusion Lyme disease model with seasonality, SIAM J. Appl. Math., 73 (2013): 2077-2099.
[69]
[ X. -Q. Zhao,
Dynamical System in Population Biology, Spring-Verlag, New York. 2003.
[70]
[ X.-Q. Zhao, Global dynamics of a reaction and diffusion model for Lyme disease, J. Math. Biol., 65 (2012): 787-808.
[71]
[ X.-Q. Zhao, Basic reproduction ratios for periodic compartmental models with time delay, J. Dyman. Differential Equations, 29 (2017): 67-82.
Lin Zhao, Zhi-Cheng Wang, Liang Zhang. Threshold dynamics of a time periodic and two–group epidemic model with distributed delay[J]. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1535-1563. doi: 10.3934/mbe.2017080
Lin Zhao, Zhi-Cheng Wang, Liang Zhang. Threshold dynamics of a time periodic and two–group epidemic model with distributed delay[J]. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1535-1563. doi: 10.3934/mbe.2017080
DownLoad:
