A systems biology strategy to identify molecular mechanisms of action and protein indicators of traumatic brain injury.

The multifactorial nature of traumatic brain injury (TBI), especially the complex secondary tissue injury involving intertwined networks of molecular pathways that mediate cellular behavior, has confounded attempts to elucidate the pathology underlying the progression of TBI. Here, systems biology strategies are exploited to identify novel molecular mechanisms and protein indicators of brain injury. To this end, we performed a meta-analysis of four distinct high-throughput gene expression studies involving different animal models of TBI. By using canonical pathways and a large human protein-interaction network as a scaffold, we separately overlaid the gene expression data from each study to identify molecular signatures that were conserved across the different studies. At 24 hr after injury, the significantly activated molecular signatures were nonspecific to TBI, whereas the significantly suppressed molecular signatures were specific to the nervous system. In particular, we identified a suppressed subnetwork consisting of 58 highly interacting, coregulated proteins associated with synaptic function. We selected three proteins from this subnetwork, postsynaptic density protein 95, nitric oxide synthase 1, and disrupted in schizophrenia 1, and hypothesized that their abundance would be significantly reduced after TBI. In a penetrating ballistic-like brain injury rat model of severe TBI, Western blot analysis confirmed our hypothesis. In addition, our analysis recovered 12 previously identified protein biomarkers of TBI. The results suggest that systems biology may provide an efficient, high-yield approach to generate testable hypotheses that can be experimentally validated to identify novel mechanisms of action and molecular indicators of TBI.

Combat injuries from the military conflicts in Iraq and Afghanistan as well as the suicides of recently retired U.S. football players have brought to light the short- and long-term consequences of traumatic brain injury (TBI; DeKosky et al., 2010; MacGregor et al., 2010). The primary insult to the head leads to secondary tissue injury, which is manifested by both immediate and delayed neurological deficits and disease (DeKosky et al., 2010; Baugh et al., 2012; Feala et al., 2013). The secondary injury process involves an intertwined cascade of evolving biochemical molecular interactions that mediate neuronal damage over hours to months after the initial trauma (Ulitsky and Shamir, 2007; Greve and Zink, 2009; Feala et al., 2013). To date, much remains unknown about how these molecular responses to injury and their poorly understood pathways are linked to clinical outcomes.

Systems biology provides an opportunity to help shed light on the extremely complex, multifactorial nature of the TBI secondary injury. Unprecedented advancements in the ability to generate high-throughput, large-scale genomic and proteomic data have led to the abstraction, construction, and graphic representation of biological networks and the ability to integrate complementary, disparate data sets in a systems-level analysis. To this end, systems biology allows for the holistic and systematic analysis of experiment-specific, high-throughput genomic and proteomic data within the context of condition-agnostic, canonical biological networks (Ideker et al., 2001). Ultimately, systems biology could be used to generate experimentally testable hypotheses.

Shojo et al. (2010) recently integrated gene expression microarray data from a fluid percussion injury (FPI) model of TBI in rats with molecular pathways to generate systems-level hypotheses and to suggest causal temporal relationships between inflammatory and apoptotic systems during the acute phase of TBI (<6 hr). These findings reaffirmed the involvement of inflammatory and survival signaling pathways that were independently observed by a different group using a systems biology approach (Kobeissy et al., 2008). The authors also hypothesized pathway links between TBI and synaptic plasticity. In another study, the function of proteins in molecular pathways was exploited to rank order and down-select potential TBI biomarkers from a list of candidates (Mondello et al., 2011).

Recently, our group reviewed strategies and potential opportunities for using systems biology to gain insights into the underlying molecular mechanisms of the TBI response and to identify novel protein indicators of brain injury (Feala et al., 2013). Through an illustrative set of 32 proteins, we showed how to use biological networks as a scaffold and overlay protein data to achieve these goals. Here, we built on these concepts and performed a meta-analysis of four large-scale gene expression data sets from distinct murine models of TBI and integrated them with biological networks to identify molecular mechanisms of action and novel protein indicators of TBI. We observed that, whereas the significantly activated biological networks were primarily associated with the immune system, the significantly suppressed networks were specific to the nervous system, particularly synaptic function. Accordingly, we hypothesized that constituent proteins from the suppressed networks were downregulated and TBI specific. Experimental testing of three such proteins by using a validated penetrating ballistic-like brain injury (PBBI) rat model of severe TBI (Williams et al., 2005) supported our hypothesis and systems biology strategy.

The strategy is based on the overarching hypothesis that certain robust molecular signatures of TBI are conserved in a meta-analysis of distinct studies involving different animal types, injury models, severity levels, and brain tissues. Implicit in this overarching hypothesis are two component hypotheses. First, a systems approach integrating gene expression data with biological network information is more likely to reveal common response mechanisms of TBI that are conserved across studies than the sole analysis of individual genes. Second, a systems approach, focused on coregulated and interconnected proteins, is effective in recovering known and identifying novel molecular mechanisms of action of TBI.

MATERIALS AND METHODS

Data Collection and Processing

In our meta-analysis, we used four publicly available gene expression microarray data sets from murine studies of moderate and severe TBI 24 hr after injury (Table1; Matzilevich et al., 2002; Natale et al., 2003; Babikian et al., 2010). These studies included two different experimental models of TBI (controlled cortical impact [CCI] and FPI), two animal types (mouse and rat), and two tissue types (cortex and hippocampus). Each study consisted of two cohorts of animals, TBI-induced animals and sham animals (i.e., controls). We downloaded the raw gene expression microarray data sets from the Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/) in September, 2011, and concurrently normalized the four data sets by using the robust multiarray average method implemented in the Bioconductor R-language suite of bioinformatics tools (Irizarry et al., 2003). This uniform standardization across the studies ensured an unbiased analysis.

Table 1

Summary of the Four TBI Gene Expression Data Sets

Data set name *	M-CCI	R-FPIm	R-FPIs	R-CCI
Reference	Natale et al. (2003)	Natale et al. (2003)	Babikian et al. (2010)	Matzilevich et al. (2002)
Model organism†	Mouse (C57BL/6)	Rat (Sprague-Dawley)	Rat (Sprague-Dawley)	Rat (Long-Evans)
TBI model	CCI	FPI	FPI	CCI
TBI severity	Moderate	Moderate	Severe	Severe
Tissue type	Cortex	Cortex	Cortex	Hippocampus
Time (hr after injury)	4, 8, 24, 72	0.5, 4, 8, 24, 72, 504	0.5, 4, 24	3, 24
Microarray platform	74Av2	U34A	U34A	U34A
Number of probes	12,488	8,799	8,799	8,799
Number of genes	7,934	4,554	4,554	4,554
Upregulated‡	667 (8.4%)	397 (8.7%)	338 (7.4%)	314 (6.9%)
Downregulated	835 (10.5%)	388 (8.5%)	372 (8.2%)	296 (6.5%)

M-CCI, mouse controlled cortical impact (CCI) model; R-FPIm, rat fluid percussion injury (FPI) model with moderate injury; R-FPIs, rat FPI model with severe injury; R-CCI, rat CCI model.

All model organisms are male.

Up-/downregulated genes were determined by using the RankProd method with a P < 0.05 cutoff (Hong et al., 2006).

Systems Biology Strategy

Figure 1 illustrates the systems biology strategy. We started by performing computational analyses to generate testable hypotheses regarding molecular mechanisms of action and protein indicators of TBI, which were then experimentally tested in the laboratory. In the computational analyses, we pursued two parallel approaches for each of the four microarray data sets. In the first approach, we separately integrated each data set with molecular pathways from a list of well-established canonical signaling and disease pathways to identify those specific to TBI 24 hr after injury. Similarly, in the second approach, we separately integrated each data set with a human protein–protein interaction (PPI) network to identify regions in the network (modules and subnetworks) representing molecular mechanisms of action specific to TBI. To generate hypotheses about molecular mechanisms and protein indicators of TBI, we considered the concordance (i.e., the conservation) of the results across the data sets for each approach and between the approaches.

Fig 1

Illustration of the systems biology strategy, in which we used computational analysis to generate testable hypotheses that can be experimentally validated in the laboratory. We started by performing a meta-analysis of four gene expression data sets and separately overlaying each data set onto two types of biological networks, canonical molecular pathways and PPI networks. In these analyses, we considered up- and downregulated genes separately because previous findings support the hypothesis that biological processes are characterized by interacting, coregulated proteins. Next, we identified statistically significant molecular mechanisms of action and protein indicators of TBI that were conserved across the studies and analyses. Finally, we experimentally tested protein indicators of TBI with an in vivo animal model. The right side of the figure illustrates the construction of a PPI subnetwork (shaded area) formed by three proteins (red circles, denoting module centers C1, C2, and C3 from three protein modules) and three other proteins (green circles) from these modules. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Because the gene expression data sets consisted of genes from rats and mice, whereas the biological networks (molecular pathways and PPI network) were represented by human genes and human proteins, respectively, for analyses we mapped the data to a common terminology, human genes, following the approach proposed by Zhang et al. (2011). We assumed a one-to-one mapping between a protein and its coding gene (i.e., we ignored protein isoforms and other translational modifications) and used Entrez Gene ID (http://www.ncbi.nlm.nih.gov/gene) to map genes and proteins, which henceforth are used interchangeably. To map genes from rat and mouse to human, we used orthology tools (Eppig et al., 2012; C. Yu et al., 2012).

Pathway Enrichment Analysis

Canonical molecular pathways provide “wiring diagrams” describing how gene products and other biomolecules interact, relate, and regulate each other to perform particular biological functions. Molecular pathways are often constructed by manual curation of literature data, which are then compiled into large databases. We collected 143 nonmetabolic canonical signaling and disease pathways representing 4,672 proteins from the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al., 2008; http://www.genome.jp/kegg/pathway.html), downloaded in December, 2011. KEGG, one of the largest and most widely used publicly available pathway databases, annotates pathways by using biological functional categories, such as nervous system, immune system, cell growth and death, and infectious diseases, which we used to identify molecular mechanisms specific to TBI. This was achieved by identifying pathways that were significantly regulated by TBI, i.e., pathways that were “enriched” with differentially expressed genes after brain injury, and associating their annotated biological functional category with brain injury.

To identify molecular pathways significantly regulated (activated by upregulated genes or suppressed by downregulated genes) at 24 hr after brain injury, we separately integrated the KEGG pathways with each of the four gene expression data sets by using the recently developed PathNet algorithm (Dutta et al., 2012). PathNet performs an extension of hypergeometric statistical test (Breitling et al., 2004), which identifies significantly regulated pathways by assessing whether the number of differentially expressed genes in the pathway is statistically significantly higher than would be expected by chance. However, unlike the hypergeometric test, which treats pathways as unordered collections of genes, PathNet capitalizes on the wiring pattern, or connectivity, among the genes in and between pathways to determine a pathway's significance within the context of the gene expression data.

PathNet separately rank ordered the association of the 143 signaling and disease pathways with TBI for each of the four data sets. We used the rank product method, a rank-based, nonparametric statistical test implemented in the RankProd routine of Bioconductor (Hong et al., 2006), to compute the statistical significance of the differential expression of the genes in each data set. The rank product method has been shown to produce differentially expressed gene lists that are more reproducible across laboratories and experimental models than other methods, especially when a small number of samples is available (Shi et al., 2006). We separately identified significantly activated (upregulated) pathways and significantly suppressed (downregulated) pathways by using a false discovery rate (FDR)-corrected P value of 0.05 as the statistical significance cutoff for both genes and pathways.

PPI Network Analysis

High- and low-throughput experimental assays are available to screen physical interactions among proteins within and between species (Rual et al., 2005; Konig et al., 2008). For interpretation and analysis, these interactions are often represented as PPI networks, in which network nodes correspond to proteins and edges between nodes represent the observed protein interactions. We constructed a comprehensive human PPI network (consisting of 74,376 physical PPIs among 11,789 human proteins) by aggregating experimentally obtained protein interaction data from nine publicly available databases (X. Yu et al., 2012), the Biomolecular Interaction Network Database (Bader et al., 2003), the Biological General Repository for Interaction Data Sets (Stark et al., 2006), the Database of Interacting Proteins (Salwinski et al., 2004), the Human Protein Reference Database (Peri et al., 2003), IntAct (Aranda et al., 2010), the Molecular Interaction database (Chatr-aryamontri et al., 2007), the mammalian PPI database of the Munich Information Center on Protein Sequences (Pagel et al., 2005), PDZBase (a PPI database for PDZ domains; Beuming et al., 2005), and Reactome (Vastrik et al., 2007). The databases were downloaded from their corresponding web sites in October, 2011.

To identify regions in the human PPI network that were significantly regulated (activated or suppressed) at 24 hr after brain injury, we separately integrated each of the four gene expression data sets with the PPI network. In particular, we searched for regions in the network where groups of connected proteins were coregulated, i.e., where groups of proteins tended to be either up- or downregulated together. We analyzed two types of PPI network regions, modules and subnetworks (constructed by combining a subset of proteins from the modules).

Module identification

For each protein in the PPI network, we defined a module as a set of connected proteins consisting of that protein (or module center Ci in Fig. 1) and its directly interacting protein partners (black circles connected to Ci). To identify TBI-specific, coregulated modules, we followed a statistical randomization procedure similar to that in our previous work detailed by Yu et al. (2011). Briefly, we identified TBI-specific modules that were either significantly activated (upregulated) or significantly suppressed (downregulated) by first computing the aggregate gene expression value of the constituent proteins in each module. Then, we randomly switched the gene expression level of each gene in all samples and recomputed the aggregate gene expression value for each module. We repeated this randomization 100 times and identified significantly activated modules (significantly suppressed modules). This was achieved by identifying those modules whose aggregate expression levels were larger (smaller) than 95% of their randomized aggregate expression levels.

Subnetwork identification

To identify larger regions of the PPI network from which to infer novel molecular mechanisms of action, we separately constructed activated and suppressed subnetworks from the corresponding statistically significant modules. In constructing subnetworks, we attempted to retain topologically important proteins from modules that were statistically significant in two or more of the four studies (i.e., modules whose significance was preserved across studies). A protein's topological importance considered the number of its interactions with other proteins, which we assumed to imply functional importance, and was used later in the downselection of protein indicators. Accordingly, we constructed interacting, coregulated subnetworks by retaining the module centers of the conserved modules and all other proteins in these modules that interacted with at least one other module center. Figure 1 provides an example of such a subnetwork (shaded area at right in the figure), constructed from three module centers (C1, C2, and C3) and three other proteins (green circles) that interacted with at least two module centers.

Subnetwork temporal profile

By using the different time points for which gene expression data were available (Table1), we calculated the temporal gene expression profile of a subnetwork as a function of time. At each time point, we computed the aggregate gene expression value of the constituent genes in the subnetwork and normalized the results such that zero represented no gene expression change, positive values indicated subnetwork activation, and negative values indicated subnetwork suppression.

Function enrichment analysis

Similarly to the pathway analysis, we also performed enrichment analysis (with the hypergeometric test) for the significantly regulated PPI modules and subnetworks. The analysis determined whether the number of proteins associated with a particular biological function in a module (or in a subnetwork) was significantly higher than would be expected by chance. However, because these network regions represent uncharacterized biological functions (unlike pathways), we used the web service resources of the Database for Annotation, Visualization, and Integrated Discovery (DAVID; version 6.7; http://david.abcc.ncifcrf.gov; Huang et al., 2009) to associate modules and subnetworks with biological functions in the form of disease, tissue, and gene ontology (GO) terms (Ashburner et al., 2000). We accessed DAVID in February, 2012, and used a significance cutoff of 0.05 (FDR-corrected P).

Known Biomarker Candidates and Downselection of Protein Indicators

We used a list of 32 previously identified TBI protein biomarker candidates (see Table1; Supp. Table S1 in Feala et al., 2013) to assess the ability of the significantly regulated pathways and subnetworks to recover such proteins. Some of these biomarker candidates have garnered multiple literature citations, are associated with TBI through diverse biological roles reported in clinical and laboratory studies, and have been found to be expressed in serum, cerebrospinal fluid, and brain tissue.

The molecular pathways and PPI network analyses generated distinct evidence for implicating protein indicators with brain injury 24 hr after the insult. Hence, to identify potential protein indicators of TBI, we considered a protein's 1) topological importance in a subnetwork, 2) presence in a significantly regulated pathway, 3) association with neurological diseases, 4) functional importance in the central nervous system, and 5) literature review. Subsequently, we experimentally validated a subset of the identified protein indicators in an in vivo animal model of TBI.

Experimental Animal Model Validation

Animals and surgical procedures

Male Sprague-Dawley rats weighing 250–300 g (Charles River Laboratories, Raleigh, NC) were used for these studies and housed individually under a normal 12-hr light/dark cycle (lights on at 6:00 am). For PBBI and sham surgery, animals were anesthetized with 5% isoflurane delivered in oxygen. The body temperature of each animal was maintained at 37°C with a heating blanket (Harvard Apparatus, Holliston, MA). Prior to biospecimen collection, animals were anesthetized with 70 mg/kg ketamine and 6 mg/kg xylaxine. Facilities at the Walter Reed Army Institute of Research (WRAIR) are accredited by the Association for Assessment and Accreditation of Laboratory Animal Care International. The experimental procedures were approved by the WRAIR Animal Care and Use Committee. Research was conducted in compliance with the Animal Welfare Act and other federal statutes and regulations relating to animals and experiments involving animals and adhered to principles stated in the Guide for the care and use of laboratory animals (National Research Council Publication, 2011 edition).

PBBI model

We separated the six rats into two equally sized groups and applied a PBBI (as previously described by Williams et al. [2006a,2006b]) to one group and a sham operation to the other group. Experimental PBBI has been extensively characterized and reproduces the temporary cavity left in the brain to mimic the ballistic nature of a high-velocity bullet wound. Briefly, a 10% unilateral frontal PBBI was induced in rats by stereotaxic insertion of a specially designed probe into the right hemisphere of the brain. The probe was inserted through a cranial window over the frontal cortex, and rapid inflation/deflation of the water-filled balloon was used to create a temporary cavity in the cerebrum. Sham rats received identical surgical procedures without balloon expansion. In both groups, animals were sacrificed 24 hr after surgery. A 2-mm slice (6 mm from bregma) of contralateral tissue and ipsilateral tissue was flash frozen in liquid nitrogen and stored at −80°C until use (three animals each for PBBI and sham).

Western blot analysis

For Western blotting, we thawed ipsilateral and contralateral tissues on ice and sonicated them in 1× radioimmunoprecipitation assay buffer (Millipore, Billerica, MA) supplemented with a protease/phosphatase inhibitor mix (Thermo Fisher Scientific, Rockford, IL). Lysate was then centrifuged at 10,000g for 20 min, and the supernatant was stored. We determined protein concentrations with the bicinchoninic acid assay (Thermo Fisher Scientific, Waltham, MA). After normalization, 10 µg of total protein from each sample was separated with 4–12% NuPAGE gels (Life Technologies, Grand Island, NY) and transferred to polyvinylidene difluoride membranes. We probed blots with antibodies for particular proteins (Abcam, Cambridge, MA; Cell Signaling Technology, Danvers, MA), and we performed densitometry and background subtraction with an LAS 4000 in Image QuantTL v 7.0 (GE Healthcare, Piscataway, NJ).

RESULTS

Significant Pathways Conserved Across Data sets

In separate analyses of up- and downregulated genes for each of the four gene expression data sets, we identified a total of 61 significantly activated pathways and 36 significantly suppressed pathways. Among these, 14 (23%) activated pathways and 19 (53%) suppressed pathways were conserved in three or more data sets (Fig. 2). These fractions of conserved pathways across the studies were much higher than the fractions obtained in a separate analysis in which we combined the up- and downregulated genes together (only one of 38 significant pathways [2.6%] was conserved in three data sets).

Fig 2

Significantly regulated pathways conserved in at least three of the four data sets.

Figure 2 lists the 33 conserved pathways, which show a sharp contrast in the functional categories between the activated (left column) and suppressed (right column) pathways. Ten (71%) of the fourteen activated pathways were related to immune response or infectious diseases, and seven (37%) of the 19 suppressed pathways were related to the nervous system. Our analysis captured seven (78%) of nine nervous system pathways represented in the KEGG database, with five pathways conserved in all four data sets and two pathways conserved in three data sets. Conversely, no activated pathway was related to the nervous system, and no suppressed pathway was related to immune response or infectious diseases.

Significant PPI Modules Conserved Across Data Sets

We identified a total of 226 significantly activated PPI modules, of which 27 (12%, consisting of 2,119 proteins) were conserved in two or more data sets (Fig. 3). Similarly, we identified 53 significantly suppressed PPI modules, of which six (11%, consisting of 296 proteins) were conserved in two or more data sets (Figs. 3, 4A). None of the 27 activated modules and only half of the suppressed modules were significant in the R-CCI data set. This finding was consistent with the results of the pathway analysis described above, in which the R-CCI study identified the lowest number of significant pathways (Fig. 2).

Fig 3

Significantly regulated protein–protein interaction (PPI) modules conserved in at least two of the four data sets.

Fig 4

A: The six suppressed PPI modules, each composed of a module center protein (red hexagon) and other proteins directly connected to the module center, for a total of 296 interacting proteins. B: Synaptic subnetwork extracted from the six modules in A consists of 58 proteins, including the six module center proteins and proteins that interacted with two or more module centers (green circles in A). Forty-five percent (26 of 58, corresponding to P < 2.0 × 10−20) of these proteins were associated with synapses and 28% (16 of 58, corresponding to P < 8.0 × 10−3) with schizophrenia. *Proteins selected for experimental testing.

To determine the biological functions associated with the conserved PPI modules, we first characterized the functions of the module center proteins. Among the activated modules, over 50% were associated with immune response or apoptosis, whereas all center proteins for the suppressed modules were associated with synaptic function (Fig. 3). We then extended the analysis and characterized the entire set of constituent proteins of the conserved modules (2,119 for the activated modules and 296 for the suppressed modules) in terms of four functional categories: GO biological process, GO cellular component, disease association, and tissue type. Table2 shows the five most enriched terms for each of the four categories in order of decreasing statistical significance (except for subnetwork entries). For the 2,119 proteins in the 27 activated modules, the functional characterization was nonspecific in each of the four categories, involving terms associated with a whole host of biological functions, including cell death. Only two entries were associated with the central nervous system (tissue type: Cajal-Retzius cell and fetal brain cortex), but each had a low prevalence (<6%). In sharp contrast, the 296 proteins in the six suppressed modules were both strongly associated with and prevalent in neurological processes and cellular components, neurological diseases, and brain tissues. For example, 21% of the proteins were associated with neuronal synapse and 66% of the proteins were expressed in brain tissues.

Table 2

Enrichment Analysis of the Proteins in the 27 Activated Protein–Protein Interaction (PPI) Modules and the Six Suppressed PPI Modules That Were Conserved in at Least Two of the Four Data Sets

Category	Activated modules*			Suppressed modules/subnetwork†
	Term			Term	Protein
		Protein modules‡			Modules§		Subnetwork§
		Count	Percentage		Count	Percentage	Count	Percentage
GO biological process
Translational elongation		84	4.0	Transmission of nerve impulse	46	15.9	18	31.0
Regulation of cellular protein metabolic process		207	9.8	Synaptic transmission	43	14.9	17	29.3
Translation		140	6.7	Cell–cell signaling	56	19.4	20	34.5
Regulation of programmed cell death		282	13.4	Regulation of synaptic transmission	27	9.3	15	25.9
Regulation of cell death		282	13.4	Regulation of neurological system process	28	9.7	15	25.9
GO cellular component
Cytosol		504	24.0	Synapse	60	20.8	26	44.8
Cytosolic ribosome		68	3.2	Synapse part	48	16.6	23	39.7
Cytosolic part		93	4.4	Cytoskeletal part	82	28.4	25	43.1
Ribosomal subunit		76	3.6	Cytoskeleton	100	34.6	26	44.8
Ribosome		96	4.6	Postsynaptic membrane	34	11.8	17	29.3
Disease association
Lupus erythematosus		38	1.8	Schizophrenia	23	8.0	16	28.0
Breast cancer		101	4.8	Huntington's disease	5	1.7	3††	5.2
Colorectal cancer		75	3.6	Cognitive function	6	2.1	2††	3.4
Crohn's disease		28	1.3	Schizophrenia/bipolar disorder	4††	1.4	3††	5.2
Ovarian cancer		37	1.8	Bipolar disorder	9††	3.1	2††	3.4
Tissue type
Cajal-Retzius cell		115	5.5	Brain	191	66.1	50	86.2
Fetal brain cortex		115	5.5	Platelet	29	10.0	3††	5.2
Epithelium		570	27.1	Fetal brain	31	10.7	10	17.2
B-cell lymphoma		68	3.2	Fetal brain cortex	16	5.5	3††	5.2
Platelet		175	8.3	Epithelium	78	27.0	13††	22.4

Involving 2,119 proteins (2,100 annotated in DAVID).

Involving 296 proteins (289 annotated in DAVID).

Number and fraction among the 2,100 proteins in the activated modules annotated by DAVID with a specific term. A protein may have multiple distinct annotations in a given category.

Number and fraction among the 289 proteins in the suppressed modules (or 58 proteins in the subnetwork) annotated by DAVID with a specific term. A protein may have multiple distinct annotations in a given category.

Italic indicates terms associated with the nervous system.

Not statistically significant (P > 0.05); for all other entries P < 0.05.

Synaptic Subnetwork

Starting with the 296 proteins in the six suppressed modules (Fig. 4A), we constructed a subnetwork by retaining the six protein module centers and 52 proteins that interacted with at least two module centers (Fig. 4B; Supp. Info. Tables I, II). When compared with the 296 proteins, the subnetwork comprising these 58 proteins was further enhanced with neurological-function-related and disease-related proteins. For example, 45% (26 of 58 corresponding to P < 2.0 × 10−20) of the proteins in the subnetwork were associated with synapses and 28% (16 of 58 corresponding to P < 8.0 × 10−3) were associated with schizophrenia, whereas the fractions of corresponding proteins in the six modules were 21% and 8%, respectively (Table2). In addition to the 26 proteins annotated as synapse, we found that the subnetwork also included three other proteins, disrupted in schizophrenia 1 (DISC1), synaptic Ras GTPase-activating protein 1, and discs large (Drosophila) homolog-associated protein 4 that play a significant role in synaptic function. Given the significant overrepresentation of proteins associated with synaptic function in this group of interacting, coregulated genes, we termed this suppressed network the synaptic subnetwork.

To investigate the functional characteristics of the synaptic subnetwork further, we analyzed its temporal expression changes from 30 min to 21 days after injury by using the R-FPIm data set (Table1). Figure 5 shows the temporal profile of the normalized aggregate gene expression values, indicating that the genes in the synaptic subnetwork tended to be attenuated during the 21-day period, reaching a nadir between 24 and 72 hr after injury and eventually returning to baseline at 21 days. A similar analysis with the more limited M-CCI data set (Table1) corroborated these results.

Fig 5

Temporal profile of the aggregate expression scores of the genes in the synaptic subnetwork (Fig. 4B) for data set R-FPIm in Table1. The gene expression was suppressed during the 21-day period, reaching a nadir between 24 hr and 72 hr after injury and eventually returning to baseline at 21 days.

Recovering Previously Identified and Inferring Novel TBI Protein Indicators

We also hypothesized that the systems biology strategy would be able to recover previously identified TBI biomarker candidates and infer novel protein indicators. To test the former hypothesis, we investigated the occurrence of the 32 proteins previously implicated with TBI in the significantly regulated pathways and subnetworks discussed above. Among these 32 proteins, 30 were represented in the human PPI network and 24 were annotated in the KEGG database. We located 12 (38%) of the 32 proteins in the significantly regulated pathways and subnetworks conserved across the studies. Table3 shows that we observed five of the 12 proteins in both analyses, eight in the PPI subnetworks (one suppressed and seven activated) and 10 in the pathways (four suppressed and six activated; see Fig. 2). We believe that this was not a chance observation because, for example, the probability of recovering at least eight of 32 proteins in an equally sized random PPI network is less than one in 5,000 (P < 2.0 × 10−4).

Table 3

Previously Identified TBI Biomarker Candidates and Novel Protein Indicators Present in the Significant Subnetworks and Pathways

Gene symbol	Gene name	Subnetwork*	Pathway†
Previously identified candidates
MAPT	Microtubule-associated protein tau§	↓	N/A‡
UCHL1	Ubiquitin C-terminal hydrolase	↑	N/A
MBP	Myelin basic protein	↑	N/A
HSPA4	Heat shock protein 70	↑	Antigen processing and presentation
CYCS	Cytochrome c	↑	Small-cell lung cancer, pathways in cancer, toxoplasmosis, p53 signaling pathway, legionellosis, tuberculosis, influenza A, herpes simplex infection
BCL-2	B-cell CLL/lymphoma 2	↑	Cholinergic synapse
			Toxoplasmosis, pathways in cancer, small cell lung cancer, tuberculosis
IL6	Interleukin 6	N/A	NOD-like receptor signaling pathway, malaria, pathways in cancer
APP	β-Amyloid (Aβ) protein 42	N/A	Serotonergic synapse
CASP7	Caspase-7	↑	Pertussis, legionellosis
CASP9	Caspase-9	↑	p53 Signaling pathway, toxoplasmosis, pathways in cancer, small-cell lung cancer
BDKRB1	B1 bradykinin receptor	N/A	Calcium signaling pathway
BDKRB2	B2 bradykinin receptor	N/A	Calcium signaling pathway
Novel protein indicators
DISC1	Disrupted in schizophrenia 1	↓	N/A
PSD95††	Post synaptic density 95	↓	Glutamatergic synapse
NOS1	Nitric oxide synthase 1	↓	Long-term depression

Significantly regulated subnetworks conserved in two or more databases (↑, activated subnetwork; ↓, suppressed subnetwork).

Significantly regulated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways from Figure 2 (italic indicates pathways associated with nervous system).

N/A, protein not present in significant pathway or subnetwork.

Tau protein is upregulated in the suppressed module.

Also known as discs large homolog 4 (DLG4).

We also used the systems biology strategy to generate testable hypotheses and infer novel protein indicators of TBI. In particular, we focused our analysis on the constituent proteins of the synaptic subnetwork, considering not only their topological importance (which could lead to truly novel biomarkers) but also distinct complementary evidence, such as their presence in significantly regulated pathways as well as their biological function and association with neurological processes. After down-selection, we arrived at three potential protein indicators of TBI that we hypothesized to be downregulated (Table3), postsynaptic density protein 95 (PSD95), nitric oxide synthase 1 (NOS1), and DISC1.

Hypothesis Testing in an Animal Model of TBI

We tested the hypothesized protein indicators at 24 hr after injury with a total of six rats separated into two groups, three subjected to 10% PBBI and three sham (controls), to determine whether abundance changes between the groups were suppressed as we had hypothesized. For both ipsilateral and contralateral tissues, we observed multiple immunoreactivity bands for DISC1 (Fig. 6A), corresponding to multiple DISC1 isoforms or their complex. Because of the indistinguishable neurological functions of these isoforms, their relative densities were compared as the aggregate of all detectable bands. Compared with sham, the abundance of DISC1 decreased by 47% in the ipsilateral tissues (P < 0.001 for between-group differences) but showed no significant differences in the contralateral tissues. We observed a significant decrease in the major isoform of PSD95 (95 kDa band) in ipsilateral PBBI tissues (69%, P < 0.05) and a nonsignificant increase (50%, P > 0.05) in the contralateral PBBI tissues (Fig. 6B). The protein's low-molecular-weight fragments (38 kDa and 25 kDa) were detectable only in the ipsilateral PBBI tissues. For NOS1, we observed an isoform of 160 kDa (α isoform) in both tissues for the two conditions, which showed a significant abundance decrease in the ipsilateral tissue (50%, P < 0.02) and the contralateral tissue (46%, P < 0.005). Another NOS1 isoform of 150 kDa (β isoform) was observed only in the ipsilateral PBBI tissues (Fig. 6C).

Fig 6

Western blot analyses of the three proteins hypothesized to have a reduced abundance at 24 hr after PBBI compared with sham (craniotomy surgery without PBBI); n = 3 animals for PBBI and sham unless otherwise noted. A: DISC1 in ipsilateral tissues shows reduced abundances (**P < 0.001). 1For ipsilateral and contralateral tissues, each bar in the graphs represents the sum of all bands on one gel lane for one animal. B: PSD95 in ipsilateral tissues also shows reduced abundances compared with sham (*P < 0.05). Unexplained increased abundances in contralateral tissues were not significant (P > 0.05). C: NOS1 shows reduced abundances in both ipsilateral tissues (*P < 0.02) and contralateral tissues (*P < 0.005). 2One PBBI sample and one sham sample yielded no signal.

In summary, we observed a statistically significant (P < 0.05) decrease in the abundance of each of the three proteins (DISC1, PSD95, and NOS1) in ipsilateral PBBI tissues (47%, 69%, and 50%, respectively) compared with ipsilateral sham tissues. Thus, we conclude that the in vivo experimental results support our hypothesis.

DISCUSSION

Our findings support the overarching hypothesis that certain robust molecular signatures of TBI are conserved in a systems biology meta-analysis of distinct studies involving different rodents and animal injury models, severity levels, and brain tissues. In particular, we showed that, by projecting high-throughput, brain gene expression data at 24 hr after TBI onto injury-independent biological network scaffolds (canonical molecular pathways and human PPI networks), we could delineate a subset of coregulated protein interactions and mechanisms of action associated with TBI that were preserved across studies. This preservation was also demonstrated in the animal study, which confirmed the downregulation of three hypothesized protein indicators in spite of using a different model of injury (PBBI) from those in the meta-analysis.

Conserved Pathways

In our functional analyses of the conserved molecular pathways, PPI modules, and subnetworks we found that the significantly activated (upregulated) responses were nonspecific to TBI and were associated primarily with the immune system, infectious diseases, and apoptosis (Figs. 2, 3, Table2), each playing important roles in mediating and resolving inflammatory responses. In sharp contrast, the significantly suppressed (downregulated) responses were associated with the nervous system in general and with synaptic function in particular. These findings are supported by previously observed experimental studies that showed an elevation of neuroinflammation in response to brain tissue damage (Shojo et al., 2010) and a reduction of synapse density as a consequence of the damage (Matzilevich et al., 2002; Natale et al., 2003). In particular, in a study of three TBI injury models, Risling et al. (2011) found that enrichment analysis of differentially expressed genes 24 hr after injury also indicated a downregulation of genes involved in neurogenesis and synaptic transmission.

Our pathway analysis showed that ∼70% of the activated pathways were related to immune response and infectious diseases (Fig. 2, left column). In fact, the six activated pathways conserved across all four gene expression studies contain one smaller pathway, toll-like receptor signaling, which is part of the innate immune system that recognizes damaged-associated molecular patterns and triggers the inflammatory response after TBI (Hua et al., 2011). Acute inflammatory response after TBI leading to neuroprotection and neurodegeneration is well known (Morganti-Kossmann et al., 2007; Blaylock and Maroon, 2011; Fahlenkamp et al., 2011), and it has been reported with the upregulation of a variety of proteins related to immunity and inflammation (Matzilevich et al., 2002).

Seven suppressed pathways were related to the nervous system, with five regulating different neurotransmitters (Fig. 2, right column). This suggests a common, nonspecific synaptic degeneration mechanism, perhaps resulting from the elevation of intracellular calcium ion that can occur through various processes and has been linked to TBI-induced neuronal loss and synaptic degeneration (Young, 1992; Bezprozvanny and Hiesinger, 2013). Nine of the twelve significantly suppressed pathways not specific to the nervous system contain a small calcium signaling pathway (Fig. 2, right column, regular font), which supports our conjecture.

Conserved Modules and Synaptic Subnetwork

In contrast to pathway analysis, which is inherently unable to reveal new protein interactions and mechanisms of action, our analysis of a large human PPI network (consisting of 74,376 interactions among 11,789 proteins) revealed the structure and composition of groups of coregulated interacting proteins. In agreement with the pathway analyses, we found that the conserved activated protein modules were associated with immune response or apoptosis (Fig. 3, left column) or were nonspecific (Table2, left column). In sharp contrast, we found all six suppressed modules to be highly associated with synaptic function (Fig. 3 and Table2, right columns). We noted that the R-CCI study was less well conserved in both the pathway and the PPI analyses (Figs. 2, 3, respectively). We believe that this was because the R-CCI study, unlike the other three, contained no replicate samples for each condition (Matzilevich et al., 2002), resulting in fewer significant differentially expressed genes (Table1).

Construction of the synaptic subnetwork from the six suppressed modules (Figs. 3, 4) delineated a highly interacting, coregulated group of 58 proteins in which 45% were associated with synaptic function and 86% were expressed in the brain (Table2; Supp. Info. Tables I, II). These findings support our hypothesis and previous observations (Chuang et al., 2007; Yu et al., 2011; Zhang and Ouellette, 2011) that proteins with similar abundances and functions tend to interact and that interacting proteins in a PPI network work together to produce a specific phenotype. Together with our recent observation that TBI biomarker candidates are highly connected in a human PPI (Feala et al., 2013), these findings allow us to speculate that some of these 58 proteins could serve as molecular indicators of TBI. Although the structure and composition of this synaptic subnetwork is likely to change with time after injury, cell type, and brain region, it allowed us to hypothesize that the downregulation of the constituent subnetwork proteins is correlated with the density of neuronal synapses and reflects the structural degeneration of neuronal synapses induced by TBI. In addition, analysis of its temporal profile (Fig. 5) supports the notion that such a degeneration of neuronal synapses is reversible. This plasticity is supported by the work of Scheff et al. (2005), who observed synaptogenesis after TBI, and the recent findings of Gao et al. (2011), who showed that synaptic density is significantly reduced at 72 hr after moderate TBI in mice and that the degree of degeneration diminishes over time. The modest activation of the subnetwork at 4 hr after injury (Fig. 5) possibly is due to the primary-injury-induced excessive glutamate release, which promotes the transcription of glutamate receptors (Wang et al., 2012) and the delayed synaptic degeneration (Gilman et al., 2003) manifested in the suppression of the synaptic network at 24 hr after injury.

We also hypothesized that the discovery of TBI-specific molecular mechanisms conserved across studies is more likely to be reproducible by using a systems biology approach than by the sole analysis of individual genes. To test this hypothesis, we computed the overlap of the constituent genes in the significant pathways and in the significant subnetworks across the four studies and compared those with the overlap of equally sized sets of top-ranked, differentially expressed individual genes from the gene expression studies. Figure 7 shows a consistent trend; the overlap in the genes extracted from the pathways and subnetworks was consistently larger than the overlap from the individual genes in the experimental list. These results support our hypothesis and hinge on the notion that, by projecting TBI-specific, high-throughput gene expression data onto injury-agnostic biological networks and essentially integrating complementary and diverse molecular information, we are better able to filter out some of the inherent variability in gene transcription data and the differences in the experimental setups. Recent studies involving comparisons across species and breast cancer data sets support our findings (Chuang et al., 2007; Zinman et al., 2011).

Fig 7

Genes in the significant pathways and in the significant subnetworks had greater overlap across the four data sets than equal-sized sets of top-ranked, differentially expressed individual genes in the data sets. We also computed the overlap of equal-sized randomly selected individual genes from each data set (1,000 realizations; error bar represents one standard deviation). A,B: Analysis of activated and suppressed pathways, respectively. C,D: Analysis of activated and suppressed subnetworks, respectively.

The ability to recover previously identified TBI biomarker candidates provided reassurance that the systems biology strategy was also capable of inferring novel protein indicators. From a list of 32 proteins previously implicated with TBI, our analysis recovered 12 proteins (38%; Table3), including the microtubule-associated protein tau, which has been shown to be predictive of clinical outcome and intracranial pressure after severe TBI (Zemlan et al., 2002; Liliang et al., 2010), and ubiquitin C-terminal hydrolase, which is currently undergoing clinical trials (Mondello et al., 2012).

Our investigation of the synaptic subnetwork along with distinct complementary evidence from the pathway functional enrichment analyses led to the hypothesis that three proteins (PSD95, NOS1, and DISC1; see Table3) would be downregulated 24 hr after TBI. From a topological perspective, both PSD95 and DISC1 were module center proteins (Fig. 2B). In addition, PSD95 formed a hub that interacted with as many as 39 proteins (67%) in the synaptic subnetwork. NOS1 directly interacted with PSD95 and served as a bridge between PSD95 and another module center, calcium/calmodulin-dependent protein kinase type IIα, an enzyme involved in calcium signaling, which is crucial for regulating glutamatergic synapses. From biological and functional perspectives, all three proteins were linked to schizophrenia (Fig. 4B; Supp. Info. Table I) and have long been studied for their associations with neurological diseases (Brenman et al., 1996; Takeuchi et al., 1997; Heales et al., 1999; Deckel, 2001; Law et al., 2001; Hodgkinson et al., 2004; Chubb et al., 2008; Brandon et al., 2009). However, to date, only PSD95 has been considered in TBI (Luo et al., 2011).

PSD95 is a scaffolding protein whose most important function is to organize glutamate receptor complexes at the postsynaptic membrane. It binds to and regulates the function and localization of a key glutamate receptor, N-methyl-D-aspartate receptor (NMDAR) subunits GRIN1, GRIN2A, and GRIN2B, which are involved in neuronal cell damage caused by excitotoxicity, a TBI-induced secondary injury (Luo et al., 2011). Because PSD95 acts as a mediator of postsynaptic signal transduction, time-dependent losses of PSD95 after TBI disrupt PSD95/NMDAR interactions and reduce neuronal cell death (Luo et al., 2011). NOS1 is a neuron-specific enzyme that synthesizes nitric oxide, a neurotransmitter involved in the regulation of synaptic signaling and synaptic plasticity (Brenman et al., 1996). The α isoform of NOS1 (160 kDa) contains a PDZ domain that binds to the PSD95/NMDAR complex, forming a complex required for cellular oxidative damage (Luo et al., 2011), which has been linked to multiple neurological diseases. The NOS1β isoform (150 kDa) lacks this domain, which is required for such an interaction. Therefore, our inferred NOS1 protein indicator should correspond to its α isoform. DISC1 was first studied in association with schizophrenia; however, its roles in Alzheimer's disease and regulation of postsynaptic density have emerged in recent studies (Brandon et al., 2009; Wang et al., 2011).

Western blot analysis of a PBBI rat model of severe TBI confirmed our hypothesis and showed that the abundance of each of the three proteins was significantly reduced (∼50%, P < 0.05) in ipsilateral tissues. The abundance reduction in PSD95 is supported by the work of Ansari et al. (2008a,2008b), who observed decreases in PSD95 levels in both the hippocampus and the cortex as early as 24 hr after injury. In contrast, for an FPI rat model, Wada et al. (1998) reported an increase in the abundance of NOS1 starting at 24 hr after injury, which became statistically significant after 3 days. This conflicts with the findings of Rao et al. (1999), who reported a significant increase in NOS1 as soon as 2 hr after CCI brain injury. We could not find previously reported results for DISC1 within the context of neurotrauma.

These results have two overarching implications. First, they show that our analysis was able to identify significantly regulated proteins regardless of experimental brain-injury model; none of the four gene expression data sets used in our meta-analysis was derived with the PBBI model used in the experimental testing (Table1). This further suggests that the systems biology strategy is robust and can detect molecular signatures of TBI independent of the mechanism of injury. Second, none of the protein-coding genes of PSD95, NOS1, and DISC1 was significantly differentially expressed in any of the four data sets (0.09 < P < 0.59). In fact, DISC1 was not even included in any of the data sets. This strongly suggests that, although genes might be unmeasured or nonsignificant in the original gene expression data set, their significance may emerge within the context of the connectivity information in biological networks.

Limitations

Our ability to perform a whole-genome analysis of brain injury in animal models was limited by multiple factors. The publicly available gene expression data sets cover only a partial number of genes (<8,000), brain regions, and time points after injury. In addition, the biological networks cover only a fraction of the nearly 20,000 human genes (∼5,000 in the pathways and ∼12,000 in the human protein-interaction network), and their construction is biased toward well-studied genes and interactions (X. Yu et al., 2012; Feala et al., 2013). Also, pathway analysis is inherently unable to reveal novel protein interactions, and algorithms to extract new molecular mechanisms of action from protein interaction networks are still immature. Our ad hoc approach to identifying protein indicators of TBI could also be susceptible to similar limitations in the discovery of novel proteins, primarily because of our reliance on established associations with known diseases and literature reviews (Kobeissy et al., 2006). We partially mitigated this by using the topological importance of a protein in the network and its presence in a significantly regulated pathway in our selection scheme, i.e., using information that is less susceptible to biases.

CONCLUSIONS

Our findings show that certain molecular signatures of TBI are conserved in a systems-level meta-analysis of four distinct gene expression studies involving different rat and mouse strains and different models of brain injury, severity levels, and brain tissues. The significantly activated molecular signatures conserved at 24 hr after brain injury were nonspecific to TBI. In stark contrast, the suppressed signatures were specific to the nervous system and associated primarily with synaptic function. We identified a suppressed synaptic subnetwork consisting of 58 highly interacting, coregulated proteins from which we hypothesized three novel protein indicators that were either not measured or not statistically significant in the original gene expression studies. We confirmed the hypothesis by Western blot analysis. This demonstrates that the importance of protein indicators can emerge within the context of biological network information. Taken together, our results suggest that systems biology may provide an alternative approach to generate testable hypotheses and identify novel molecular mechanisms of action and protein indicators of TBI.