Use of quasi-normal form to examine stability of tumor-free equilibrium in a mathematical model of bcg treatment of bladder cancer

1.   Introduction

Traumatic brain injury (TBI) results from external forces such as a blow to the head or an object penetrating the skull. TBI is regarded as a common cause of death and disability and a major healthcare issue worldwide ^[1,2]. It often leads to various central nervous system (CNS) symptoms and greatly increases the risk of developing long-lasting psychiatric disorders. Since TBI can be highly heterogeneous, the prevention of post-TBI complications is one of the greatest challenges for researchers and therapeutists.

Neuroinflammation is considered as an important post-TBI secondary injury mechanism, with hallmarks such as neurovascular unit (NVU) alterations ^[3]. The NVU is mainly composed of neurons, astrocytes, microglia, endothelial cells and smooth muscle cells ^[4], and it has been demonstrated that the activation of microglia and astrocytes could promote neuroinflammation via secreting more proinflammatory cytokines, such as TNF-α and IL-1β ^[1,5,6]. Certain microglia receptors, namely CX3CR1, P2Y12R, TLR4, IL1R1, RAGE and TNFR1, have been highlighted, and their associated downstream pathways include the PI3K/AKT, MAPK, cAMP, and NF-κB signaling pathway, which activate the expression of a series of inflammatory genes ^[7]. Previous studies have also suggested an association between neuroinflammatory responses and excitatory-to-inhibitory synaptic imbalance in patients and animal models of TBI. Synaptic reorganization has been observed in mice after TBI, which is possibly caused by disturbed neurotransmitter release ^[8,9]. In addition, TBI could lead to endothelial dysfunction and cardiovascular diseases, where the nitric oxide (NO) and endothelial-dependent hyperpolarization (EDH) pathways played a critical role ^[5]. Meanwhile, vascular endothelial growth factor (VEGF) has been found to participate in neuroprotection after TBI ^[10].

The advances in single-cell RNA sequencing technology have enabled the exploration of complex and rare cell populations and helped uncover regulatory mechanisms between genes and cells involved in the pathogenic processes of several diseases ^{[11,12,13,14]}. However, a detailed landscape of molecular mechanisms and pathways involved in post-TBI pathology still remains absent. In the present study, we have integrated single-cell and bulk RNA-seq data to identify TBI-related cell types and regulators, hopefully providing some novel therapeutic targets for TBI.

2.   Materials and methods

2.1.   Data acquisition

The single-cell RNA-seq data from cortex of mice was collected from Gene Expression Omnibus (GEO) with accession number GSE124952 ^[15]. The RNA-seq data of cortex tissues in mice with TBI and normal controls and the RNA-seq data of Dap12 (Tyrobp)-deficient mouse microglial cells were downloaded from GEO with accession numbers GSE79441 ^[16] and GSE9043, respectively. Data from GEO under accession numbers GSE68207 ^[17] and GSE155063 ^[18] were used for TYROBP validation. All these data were pre-normalized in GEO database, and were used for downstream analysis.

2.2.   Cell clustering analysis

The cell clustering analysis was implemented in R Seurat ^[15] package. Specifically, the count-based expression data was first imported into R using CreateSeuraObject, and the top 2000 highly variable features were selected using FindVariableFeatures with vst method. We selected different clustering resolutions at 0.1, 0.3, 0.5 and 0.7, manually compared those clustering results, and observed that the clustering with a resolution at 0.3 achieved the optimal partitions. Thus, the clusters were identified at a resolution of 0.3 using FindClusters, and T-distributed Stochastic Neighbor Embedding (t-SNE) was applied to reduce the dimensionality. The cell-type specific genes were detected using FindAllMarkers function at adjusted P-value < 0.05 and minimal percentage > 0.25. Reported marker genes of those cell clusters were collected from an earlier study ^[16]. Visualization tools for the scRNA-seq data such as t-SNE, heat-map, and vioplot were also implemented in R Seurat package.

2.3.   Gene set enrichment analysis

Two enrichment algorithms, gene set enrichment analysis (GSEA) and gene set overrepresentation enrichment analysis (ORA) were employed in this study. These two algorithms were applied to identify cell types involved in TBI. Also, GSEA was applied to identify pathways enriched in the gene sets of interest, which was implemented in R clusterProfiler ^[17] package.

2.4.   Cell-cell communication

The ligand-receptor pairs were collected from a public dataset provided by an earlier study ^[18]. The cell-type specificity of differentially expressed genes encoding ligands or receptors in TBI were then tested in identified cell populations. Subsequently, these cell-type specific ligands or receptors involved in TBI were mapped to known ligand-receptor pairs, and those successfully mapped pairs were considered as links between corresponding cell types.

2.5.   Differential expression analysis

Differential expression analysis of bulk RNA-seq was conducted using R DESeq2 package. Specifically, disease status (TBI or control) was regarded as a variable in the meta-data. Genes with an FPKM value above 1 in at least 20% of samples were retained for this analysis. The count-based gene expression was normalized using the median of ratios method. The thresholds of adjusted P-value and log2 (fold change) were selected at 0.05 and 0.5, respectively. Differential expression analysis of single-cell RNA-seq data based on Wilcoxon rank sum test was performed, which was implemented in R Seurat package with function FindAllMarkers.

3.   Results

3.1.   Characterization of cell clusters in mice cortex by five previously reported cell types

To define cell types in mouse cortex, we collected a total of 10,538 single cells in cortex without any treatments from a previous study ^[16], and conducted T-distributed Stochastic Neighbor Embedding (t-SNE) analysis to identify cell clusters. Totally, we identified 13 cell clusters and characterized cell types by known markers ^[16] of astrocytes (Astro), oligodendrocyte, newly formed oligodendrocytes (NF Oligo), oligodendrocyte precursors (Oligo), microglia (Microglia), endothelial cells (Endo), excitatory and inhibitory neurons (Excitatory and Inhibitory) (Figure 1A). Particularly, excitatory and inhibitory neurons and endothelial cells could be further divided into 5, 2 and 2 sub-clusters in this clustering (Figure 1B), and these cell populations were denoted as Excitatory-1/2/3/4/5, Inhibitory-1/2 and Endo-1/2. Specifically, marker genes including Ctla2a, Ly6c1, Pglyrp1, Cldn5 and Itm2a were shared by these two types of endothelial cells, while Myl9, Higd1b, Cald1, Rgs5 and Vtn were specifically expressed in Endo-1. Notably, Myl9 and Cald1 were involved in vascular smooth muscle contraction. Comparison of the representative genes in the excitatory neurons revealed that gene expressions in Excitatory-1, 2 and 4 shared a high similarity, and Excitatory-3 and 5 had more unique expression patterns. Specific expressions of certain marker genes further confirmed that this characterization of cell types makes biological sense (Figure 1C), including Slc17a7, Gja1, Aspa, Bmp4, C1qa, Flt1 and Gad2, which were detected respectively in excitatory neurons, astrocytes, oligodendrocyte, newly formed oligodendrocytes, microglia, endothelial cells, and inhibitory neurons.

Figure 1.  The molecular characteristics of cell types in mouse cortex. (A) The t-SNE-plot showing the distribution of samples in the different cell-types. (B) The expression patterns of top-five marker genes in the 13 cell types. (C) The vioplots for marker genes and cell types.

DownLoad: Full-Size Img PowerPoint

3.2.   Identification of cell types involved in traumatic brain injury

To identify cell types involved in traumatic brain injury (TBI), we conducted a systematic data analysis by integrating scRNA-seq and bulk RNA-seq data. The bulk RNA-seq data was collected from a previous study ^[19]. Specifically, differential expression analysis of mouse cortex RNA-seq data revealed 2977 upregulated and 1726 downregulated genes in TBI as compared with controls (TBI (n = 3) vs. control (n = 3), adjusted P-value < 0.05 and log2 (fold change) > 0.5). Gene set overrepresentation enrichment analysis (ORA) of the upregulated and downregulated genes against the cell-type specific marker genes by scRNA-seq data revealed that the upregulated genes were highly enriched in microglia and endothelial cells, while the downregulated genes were enriched in excitatory neurons (Figure 2A, FDR < 0.05). Consistently, gene set enrichment analysis (GSEA) further confirmed that these cell types were either highly activated or inactivated in TBI (Figures 2B-2H, FDR < 0.05). These results indicated that the systematic analysis of scRNA-seq and bulk RNA-seq data could accurately detect various cell types involved in TBI.

Figure 2.  The activities of the cell types in TBI. (A) The cell types significantly altered in TBI by overrepresentation enrichment analysis (ORA). The rows represent the cell types, and the columns represent the upregulated and downregulated gene sets. The diameter of the circles represents the ratio of the number of upregulated or downregulated genes to that of cell type-specific genes. The expression levels of the cell type-specific genes were showed in (B–H). The genes ranked on the left and right represent the upregulation and downregulation in TBI. The enrichment of differentially expressed genes by TBI in the microglia (B), endothelial cells (C, D) and excitatory neurons (E–H) was tested by gene set enrichment analysis (GSEA).

DownLoad: Full-Size Img PowerPoint

3.3.   Cell-cell communications between activated cell types in TBI

To predict cell-cell communications between those activated cell types, we linked those cell populations by ligand-receptor interactions. We first identified those upregulated genes encoding ligands or receptors in each cell type, and then mapped the receptor-ligand pairs according to the curated ligand-receptor interactions from an earlier study ^[18]. As shown in Figure 3A, autocrine ligand-receptor interaction was observed in microglia and endothelial cell (Endo-1) via Icam1-Itgam and Col4a1/2-Cd93 pairs. Moreover, microglia could also communicate with endothelial cells via Icam1-Il2rg and C1qa-Cd93.

Figure 3.  The cell-cell communication network and downstream pathways activated in TBI. (A) The cell-cell communication network between the activated cell types in TBI. The cells expressing ligands and receptors were displayed in the first and the last columns. The ligands/receptors secreted by the cells expressing ligands or receptors are directly connected with their cells. The connected ligands and receptors mean that they had the potential to interact with each other. (B) The autocrine ligand-receptor of microglia and intercellular ligand-receptor communication between microglia and endothelial cells. The gene symbols on the bottom represent the genes involved in the two pathways, TYROBP causal network and integrin signaling, respectively.

DownLoad: Full-Size Img PowerPoint

Furthermore, functional enrichment analysis of the upregulated marker genes in microglia revealed that autocrine ligand-receptor interaction in microglia via Icam1-Itgam might result in activation of TYROBP causal network (Figure 3B, FDR < 0.05, C1qc, Tyrobp, Plek, Apbb1ip, Itgam, Hcls1, Cd84, Ncf2 and Pycard). The intercellular communication between microglia and endothelial cells via C1qa-Cd93 might activate the integrin signaling (Figure 3B, FDR < 0.05, Fn1, Vwf and Tln1).

3.4.   Tyrobp may be a candidate therapeutic target for TBI

As chronically activated microglia could release pro-inflammatory molecules, resulting in further tissue damage and contributing potentially to neurodegeneration in TBI, we investigated whether the microglia activation could be attenuated by inhibiting TYROBP causal network, especially through targeting the key regulator, Tyrobp, in this network. Fortunately, we found a gene expression dataset of three wild type and three Tyrobp/Dap12-deficient mouse microglia. As expected, downregulated genes in Tyrobp/Dap12-deficient microglia were highly enriched in microglia activation (Figure 4A). Specifically, 42 genes involved in microglia activation were upregulated in TBI, but were downregulated in Tyrobp-deficient microglia as compared with the controls (Figure 4B). To further confirm the upregulation of Tyrobp in TBI, we collected two gene expression datasets from TBI-affected and control animal models. Consistently, Tyrobp was confirmed to be upregulated in TBI-affected samples of the two datasets (Figure S1). These results indicated that Tyrobp deficiency could significantly reduce the microglia activity, and Tyrobp might be a candidate therapeutic target for TBI so as to inhibit pro-inflammatory molecules and reduce tissue damage.

Figure 4.  The microglia activation and signature genes in TBI. (A) The enrichment degree of genes involved in microglia activation in Tyrobp-KO mice. (B) The signature genes of microglia activation upregulated in microglia of TBI and downregulated in Tyrobp-deficient microglia.

DownLoad: Full-Size Img PowerPoint

3.5.   Cell-cell communications between inactivated cell types in TBI

As the excitatory neurons were inactivated in TBI, we then investigated the cell-cell communications between those inactivated cell types. In accordance with the communications between activated cell types, the autocrine ligand-receptor pairs were also observed in type-3 excitatory neurons (Excitatory-3, Figure 5A, FDR < 0.05). Furthermore, functional enrichment analysis of those downregulated marker genes in excitatory neurons revealed that terms including PDZs pathway, neurexins and neuroligins, disruption of postsynaptic signaling by CNV, protein-protein interaction at synapses, and neuron system were significantly downregulated in Excitatory-3 (Figure 5B). Particularly, the receptor Nrxn3 was involved in these pathways (Figure 5C), suggesting that downregulation of these ligands and receptors might cause dysfunction of nervous system-related pathways, thereby impairing the normal brain function after TBI.

Figure 5.  The cell-cell communication network and downstream pathways inactivated in TBI. (A) The cell-cell communication network between the inactivated cell types in TBI. (B) The pathways enriched by the downregulated signature genes of excitatory neurons in TBI. (C) The representative marker genes in excitatory neurons of TBI.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

In this study, we aimed to identify key cell types and regulators involved in traumatic brain injury. Specifically, we have successfully identified and characterized 13 cell populations, including astrocytes, oligodendrocyte, newly formed oligodendrocytes, microglia, two types of endothelial cells, five types of excitatory and two types of inhibitory neurons (Figure 1A). The cell-type specific expression of marker genes, including Slc17a7, Gja1, Aspa, Bmp4, C1qa, Flt1 and Gad2, which were found respectively in excitatory neurons ^[20], astrocytes ^[21], oligodendrocyte ^[22], newly formed oligodendrocytes ^[23], microglia ^[24], endothelial cells ^[25], and inhibitory neurons ^[26], confirmed the validity of proposed characterization of cell types.

Differential expression analysis and GSEA revealed that marker genes of microglia and endothelial cells were among the upregulated genes, while those of excitatory neurons were among the downregulated genes in TBI. Microglia activation has been observed in both animal models ^[27] and TBI patients ^[28], and can persist for several years ^[29]. When TBI occurs, the microglia can migrate into the sites of injury, and establish a protective environment mitigating the damage of TBI ^[30]. However, long-term activation of microglia can release pro-inflammatory cytokines and cause tissue damage ^[31]. These results indicated that microglia activation in TBI may mitigate the damage of TBI in the short term, but may cause tissue damage by releasing pro-inflammatory cytokines in the long term.

The cell-cell communication analysis revealed that microglia and endothelial cell (Endo-1) might interact via Icam1-Il2rg and C1qa-Cd93, and microglia might also communicate with each other via Icam1-Itagm. Particularly, autocrine ligand-receptor interaction in microglia via Icam1-Itgam might result in activation of TYROBP causal network (Figure 3B, FDR < 0.05, C1qc, Tyrobp, Plek, Apbb1ip, Itgam, Hcls1, Cd84, Ncf2 and Pycard). The receptor ITGAM is required, along with TYROBP, in microglia to control the pro-inflammatory cytokines ^[32]. Furthermore, as a ligand of C1q, C1QA has been widely reported to be bounded to CD93 in endothelial cells ^[33,34]. CD93 could promote β1 integrin activation during tumor angiogenesis ^[35], indicating that the ligand-receptor interaction between C1QA-CD93 could promote integrin signaling activation in endothelial cells. As chronically activated microglia could release pro-inflammatory molecules ^[31], resulting in further tissue damage and potentially contributing to neurodegeneration in TBI, we further investigated those downregulated genes in Tyrobp/Dap12-deficient microglia, and found that those genes were highly enriched in microglia activation, suggesting the upregulation of Tyrobp and TYROBP causal network in microglia might be a candidate therapeutic target in TBI ^[36,37]. In contrast, the excitatory neurons might help maintain the normal brain function ^[38], and their inactivation might cause dysfunction of nervous system in TBI patients.

5.   Conclusions

In conclusion, the present study has identified key cell types such as microglia, endothelial cells and excitatory neurons, and revealed key regulator in TBI such as TYROBP, C1QA and CD93, which shall improve our understanding of the pathogenesis of TBI.

Acknowledgments

This work was supported by The National Natural Science Foundation of China (No. 81871458), Talents Training Program of Pudong Hospital affiliated to Fudan University (Project no. PX202003), and the main subject of the discipline construction of the Pudong New Area Health and Family Planning Commission (PWZzk2017-16).

Conflict of interest

The authors declare that they have no conflicts of interest.

Supplementary

Figure S1.  The upregulation of Tyrobp in TBI-affected samples confirmed by two independent datasets. (A and B) display the expression levels of Tyrobp in TBI and controls. The red and blue boxes represent the TBI and controls, respectively.

DownLoad: Full-Size Img PowerPoint