Temporal Genetic Modifications after Controlled Cortical Impact—Understanding Traumatic Brain Injury through a Systematic Network Approach

Traumatic brain injury (TBI) is a primary injury caused by external physical force and also a secondary injury caused by biological processes such as metabolic, cellular, and other molecular events that eventually lead to brain cell death, tissue and nerve damage, and atrophy. It is a common disease process (as opposed to an event) that causes disabilities and high death rates. In order to treat all the repercussions of this injury, treatment becomes increasingly complex and difficult throughout the evolution of a TBI. Using high-throughput microarray data, we developed a systems biology approach to explore potential molecular mechanisms at four time points post-TBI (4, 8, 24, and 72 h), using a controlled cortical impact (CCI) model. We identified 27, 50, 48, and 59 significant proteins as network biomarkers at these four time points, respectively. We present their network structures to illustrate the protein–protein interactions (PPIs). We also identified UBC (Ubiquitin C), SUMO1, CDKN1A (cyclindependent kinase inhibitor 1A), and MYC as the core network biomarkers at the four time points, respectively. Using the functional analytical tool MetaCore™, we explored regulatory mechanisms and biological processes and conducted a statistical analysis of the four networks. The analytical results support some recent findings regarding TBI and provide additional guidance and directions for future research.


Introduction
Traumatic brain injury (TBI) is not only a primary injury, caused by external physical forces, but also a secondary injury caused by a series of biological processes. These include metabolic, cellular, and other molecular events that eventually lead to the death of brain cells, tissue and nerve damage, and atrophy [1]. Due to the importance of the brain, damage to this organ usually causes particularly complex symptoms relative to other injuries. For example, if the blood-brain barrier (BBB), a selectively permeable barrier to the transport of crucial elements for the brain, is damaged by an external force, several complications can occur due to disrupted homeostasis of the brain [2]. In order to take all aspects of the consequences of a TBI into account, and as knowledge on TBIs develops, treatment has become extremely complex and increasingly difficult [3]. Immediately after a TBI, primary injuries resulting from the initial trauma cause rupture of blood vessels, hemorrhaging, structural deformation of tissues, and massive death of neurons and other cells at injured sites [4]. Treatment of the primary injuries of a TBI is similar to those of other traumas. As a result of advances in the treatment of emergency traumas, secondary injuries have become the main cause of morbidity and mortality in TBI patients. Several studies have therefore concentrated on the treatment of secondary injuries resulting from TBIs in an effort to reduce the extent and effect of these injuries.
Secondary injuries are influenced by the vascular perturbations, cerebral metabolic dysfunction, and inadequate cerebral oxygenation related to primary injuries [5]. Vascular perturbations can affect the permeability of the BBB, pinocytotic activity of endothelial cells, and the redistribution of ions and neurotransmitters [6,7]. Increased permeability of the BBB leads to edema and subsequent brain swelling, and can increase the risk of brain infections [8]. Intracellular accumulation of potassium and calcium ions can affect the functioning of neurons in information transduction and of mitochondria in metabolism [9]. Elevated glutamate can be observed as early as 5 min after experimental trauma, which may cause excitotoxicity after a TBI [10] and can also affect mitochondrial functioning through oxidative stress [6]. Increased oxidative stress following a TBI is directly related to the pathogenesis of TBIs. Several oxidative markers, including reduced glutathione (GSH), the GSH/oxidized glutathione (GSSG) ratio, glutathione peroxidase (GPx), glutathione reductase (GR), glutathione-S-transferase (GST), glucose-6-phosphate dehydrogenase (G-6PD), superoxide dismutase (SOD), and catalase (CAT) are observed after a TBI. Hence, development of antioxidant strategies is of primary interest in ongoing efforts to optimize brain injury treatment [11]. Recently, inflammation was also recognized as a critical element in recovering from a TBI provided that TBI treatment is an option [12]. Inflammation is a general response to external insults and is also modulated by oxidative stress, mitochondrion-mediated metabolism, and cellular states [5]. An interwoven network of a wide variety of functions complicates TBIs and impedes the development of efficient treatments. Therefore, a systematic perspective of the evolution of a TBI is needed to provide a global view of molecular interactions.
In this study, we used a systems biology approach to elucidate protein-protein interactions (PPIs) very soon (4 h) and longer (8,24, and 72 h) after a TBI [13]. A systems biology approach can provide a systematic view on the interwoven network of functions and interactions between the proteins involved in these complex networks, which are regarded as systems with sub-units (e.g., proteins) connected as a whole. Many PPIs can form PPI networks (PPINs), in which proteins are nodes and their interactions are the edges of the network. We also used a mathematical model with the microarray data to generate quantitative descriptions of the molecular interactions. By comparing injured and control networks, we were able to select some specific interactions and proteins as potential therapeutic or monitoring targets, according to their differential interactions. Connections between these targets and the physiological phenomena of a TBI provide insights into early interventions for TBIs.

Results and Discussion
We built the network biomarkers for the four time points post-TBI. We then used the commercial pathway analysis software MetaCore™ and the free network ontology analysis (NOA) web server to do a functional pathway analysis to reveal the underlying molecular mechanisms of TBI. The main results of this research are the network biomarkers generated by our model, i.e., the proteins with the highest TBI relevance values (TRVs), and the network structure. However, MetaCore™ and NOA enabled us to add deeper medical and biology significance to our results and make them useful for identifying novel strategies for therapy or recovery processes. The powerful MetaCore™ software uses data mining and statistical methods to select valuable information from the biological and medical literature. It provides a more general perspective and enhanced interpretation of our research. The original results from our primary analysis thus form the core findings of this study, and the additional results generated by MetaCore™ and NOA can be viewed as supplementary information that is of wider benefit.

Evolution of Network Biomarkers at the Four Post-Traumatic Brain Injury (TBI) Time Points
We built a Differential PPIN (DPPIN) for each of the four post-TBI time points (4,8,24, and 72 h) ( Figure 1) and calculated the TRV of each protein in the four networks ( Table 1). The network diagrams ( Figure 1) contain more information than just the TRVs, such as that encoded by the edges and nodes. Node size represents the TRV of each protein, and edge width the interaction ability between the two proteins. Red and blue edges respectively indicate positive and negative values of d ij in Equation (7). We identified statistically significant network marker proteins for the four post-TBI stages by screening the TRV p-values. As in our previous study on stroke, we wanted to reveal the repair mechanisms that operate at these four post-TBI time points. After fold change screening (FC (Fold Change) > 1.5), we identified 27, 50, 48, and 59 significant proteins at 4, 8, 24, and 72 h post-TBI, respectively. We only list the top 20 for each time point in Table 1, but provide the full lists in a supplementary table. Their corresponding TRVs were in the ranges 4.5-64.5, 4.9-17, 5.3-30.4, and 5.1-34.2, respectively. We used these significant proteins and their PPIs to construct network markers for the four post-TBI time points. The protein relevance values for TBI were much smaller than those we calculated for carcinogenesis in our previous study on cancer [9,10], and the cancer networks were much more complex than the TBI networks. However, the TBI values and network structures were similar to those we identified for stroke. We do not discuss the UBC (Ubiquitin C), protein in this paper because it is a complex issue. UBC is a housekeeping gene for many different diseases, such as cancers, stroke, and TBI. Recent findings also demonstrate that Ubiquitin Carboxyl-Terminal Esterase L1 (Ubiquitin Thiolesterase, UCHL1) is a promising biomarker that physically interacts with UBC [14,15]. Our demonstration of UBC's role in the first 72 h after TBI supports the significance of UCHL1 in TBI. We also identified four significant proteins of core network biomarkers (intersections) for the four post-TBI time points: UBC, SUMO1, CDKN1A, and MYC (Table 2). Core and specific network biomarkers are important when discussing evolutionary processes. The discovery of these significant proteins, which are consistent throughout the four post-TBI time points, is especially interesting as Natale et al. [16] show MYC was the only one of these proteins identified as common to the different TBI models using controlled cortical impact (CCI) and fluid percussion injury (FPI). To comprehensively discuss the evolutionary behavior of TBI, we also summarized the core network (intersections) of each adjacent time point, i.e., 4-8, 8-24, and 24-72 h. We also identified four significant proteins of core network biomarkers (intersections) for the four post-TBI time points: UBC, SUMO1, CDKN1A, and MYC (Table 2). Core and specific network biomarkers are important when discussing evolutionary processes. The discovery of these significant proteins, which are consistent throughout the four post-TBI time points, is especially interesting as Natale et al. [16] show MYC was the only one of these proteins identified as common to the different TBI models using controlled cortical impact (CCI) and fluid percussion injury (FPI). To comprehensively discuss the evolutionary behavior of TBI, we also summarized the core network  Red edge means the difference in the interaction activity ≥ mean + STD; Yellow node means the protein with significant CRV. Of the biomarkers identified in our analysis, UBC has also been identified in other systems biology reviews of TBI datasets. As Feala [17] mentions, "Animal experiments do not always reproduce the same results across studies. This is primarily because of variations in animal species, injury type and severity, time course of collection, and sampled tissue". For example, MAPT (Microtubule Associated Protein Tau), also known as Tau (Tau proteins are proteins that stabilize microtubules), appears to be a promising biomarker, yet it does not appear to increase in mTBI  Of the biomarkers identified in our analysis, UBC has also been identified in other systems biology reviews of TBI datasets. As Feala [17] mentions, "Animal experiments do not always reproduce the same results across studies. This is primarily because of variations in animal species, injury type and severity, time course of collection, and sampled tissue". For example, MAPT (Microtubule Associated Protein Tau), also known as Tau (Tau proteins are proteins that stabilize microtubules), appears to be a promising biomarker, yet it does not appear to increase in mTBI (mild traumatic brain injury) scenarios. The change in genetic profiles across different time points was especially problematic for us because a practical biomarker should ideally be present and examinable at different time points after the injury. Of the four proteins common among the different time points, only UBC was mentioned in Feala's article [17]. We feel that other proteins that are persistent throughout the course of sampling should also receive additional attention.
The unique genetic profile for UBC, SUMO1, and CDKN1A at the four time points for CCI suggests the need to further assess current biological modeling designs. We have previously found that even the sham control groups require more precise treatment in CCI experimentation [17] and have proposed that unicortical drilling should be performed on the control group. Since improved biological modeling requires additional funding and repeating previous experiments, we have used a systems biology approach to re-analyze the controversial penumbra region. Given that the penumbra region contains salvageable brain tissue and that it lies between the heavily-injured region and the nearly-normal region, we feel that appropriate analysis of this region should disclose valuable information and targets for treatment.
Our analyses revealed significant information regarding the recovery and damage processes in the post-TBI stages. This information provides clues for selecting novel drug targets for therapeutic and recovery processes. Because of article length restrictions, we present these results in a supplementary table. Figure 1 and Table 1 show the difference between the network structures of the four time points. At 4 h post-TBI, the UBC node dominated the network (TRV = 64.8), while the TRVs of the other nodes were all <15. We discuss the next two most important proteins, heat shock protein HSP90AA1 (Heat shock protein 90) and ITGA4 (integrin α-4 precursor gene) (TRV = 12.4 and 11.9, respectively), below. At 8 h, there were no TRVs greater than 20. However, ELAVL1 (embryonic lethal abnormal vision-like protein 1) (TRV = 17.0) appeared to be a significant protein of neuron-related diseases. At 24 h, CDK2 (cyclin-dependent kinase-2), FN1, and epithelial growth factor receptor (EGFR) (TRV = 30.4, 25, and 20, respectively) dominated the network. We discuss these significant nodes below. Finally, at 72 h, APP (amyloid precursors protein) and ELAVL1 (TRV = 34 and 31.9, respectively) dominated the network. The color and width of the edges in Figure 1 reveal some information about the regulatory mechanisms of the proteins, but the model we used in this research focused on PPIs, not genetic regulatory networks (GRNs). We think that the significant nodes we identified are more like "hubs" than drivers. Our team has also developed GRN models [18], but due to limited data types it was not appropriate to apply a GRN model to this dataset.

Network Structure Interpretation of the Four Post-TBI Time Points
We discuss the five aforementioned proteins here.
(v) EGFR: This is a well-known cancer oncogene and the vascular EGFR (VEGFR) receptor is always reported in brain injuries [22]. It is a novel clue for us to identify the relationship between EGFR and TBI.

Comparison with Our Previous Results for Stroke
We compare our results with our recently conducted work on human stroke [23]. We identified five significant proteins common to both stroke and TBI: UBC, APP, NEDD8, PAN2, and EVAVL1.
APP activates voltage-dependent calcium channels and may induce neuronal apoptosis. Its protein product binds growth factor receptor and plays a role in the regulation of peptidase activator and acetylcholine receptor activities. It is involved in the positive regulation of peptidase activity, locomotor behavior, and axon cargo transport. It also participates in the glypican signaling pathway and Alzheimer disease (AD) pathway [24,25]. It is involved in many cellular behaviors and thus forms one of the key hubs in TBI.
NEDD8 (neural precursor cell expressed) encodes a protein that exhibits ubiquitin protein ligase binding. It is involved in protein neddylation. It participates in the p53 signaling pathway, the neddylation pathway, and the cullin-dependent proteasome degradation pathway. It is also associated with Parkinson's disease [25,26].
There are few reports that PAN2 is directly related to TBI or stroke, so this could be a novel target for therapy. The ELAVL1/Hu family of RBPs (RIM binding proteins) plays a key role in neuroscience. Skliris et al. [27] discuss how neuroprotective behavior requires the functions of the RNA-binding protein, HuR, and give a full description of the mechanisms of ELAVL1/HuR. Usually, three members of this family (HuB/HEL-N1/ELAVL2, HuC/ELAVL3, and HuD/ELAVL4) are expressed by neurons, but ELAVL1/HuR is always expressed in both neuronal and non-neuronal tissues. These findings concur with descriptions of the involvement of synaptic dysfunction and cytoarchitectural degradation in stroke [28].

MetaCore™ Results
We used the MetaCore™ Analyze Networks (AN) algorithm to analyze the network biomarkers that were unique to or common to all or some of the networks we identified for the four time points post-TBI (Table 3). The detailed algorithm is given in Supplementary Material S0, showing that completely different cellular behaviors are operating at the four time points. Briefly, most biomarkers at 4 h post-TBI were related to cell activation and signaling behaviors. In contrast, most at 8 h post-TBI were related to immune response. At 24 h post-TBI, most were related to DNA behavior, whereas in the final 72 h post-TBI stage, most were related to regulation and apoptotic behaviors. The evolution of cellular behavior at these four time points post-TBI is thus clear.
A canonical pathway in MetaCore™ represents a set of consecutive signals, or metabolic transformations, confirmed as a whole by experimental data or by inferred relationships. They are linear sets of carefully defined steps that form a map; in MetaCore™, complete biochemical pathways or signaling cascades in the commonly accepted sense correspond to maps. Over 70,000 pathways are stored in MetaBase™, which are mainly used for network generation and visualized as canonical pathway maps. The AN algorithm creates a large network and divides it into smaller sub-networks that can be separately built. Networks are more interactive than pathway maps and can sometimes be more complex. Table 3. Analysis of networks using MetaCore™. The gene content of the uploaded files was used as the input list for generating biological networks using the Analyze Networks (AN) algorithm with default settings. This is a variant of the shortest paths algorithm, with the following main parameters: (1) relative enrichment with the uploaded data and (2) relative saturation of networks with canonical pathways. These networks are built on the fly and are unique to the uploaded data. In this workflow, the networks are prioritized based on the number of fragments of canonical pathways in the network. S: size; T: target; P: pathways; G: gScore.   i.e., four). They overlap where they have network objects in common. The different elements of the resulting Venn diagram comprise of the following: (1) a "common" area in the middle, including the network objects found in all the experiments (in this case, at all four time points); (2) an area "unique" to each experiment, consisting of objects unique to that experiment; and (3) one or more sets of objects common to two or more experiments, but not all (MetaCore™). Figure 2a gives the top 10 pathway maps. Figure 2b gives the top 10 process networks. They are highly related to the cell cycle, DNA damage, and signal transduction. These pathways and process networks may provide clues for new uses of conventional drugs. From the literature review, we also found that it was important to pay attention to the relationship between the cell cycle and TBI [29]. Simone et al. discuss the fact that cell cycle inhibition provides neuroprotection and reduces glial proliferation and scar formation after traumatic brain injury. This article has been cited more than 250 times, and therefore the discussion of the relationship of the cell cycle and TBI is very important. Figure 2c gives results for related diseases. While one can make comparisons of common hidden embryological and molecular mechanisms of these diseases, the correlation and embryological relationship amongst TBI and these diseases, however, would require additional investigation and validation.   (a) Pathway maps. Canonical pathway maps represent a set of signaling and metabolic maps that comprehensively cover human biological functioning. All maps are created by Thomson Reuters scientists using a high quality manual curation process based on published peer reviewed literature. Experimental data are visualized on the maps as blue (for downregulation) and red (upregulation) histograms. The height of the histogram corresponds to the relative expression value for a particular gene/protein; (b) Process networks. The content of these cellular and molecular processes is defined and annotated by Thomson Reuters scientists. Each process represents a pre-set network of protein interactions characteristic of that process; (c) Related diseases (determined based on biomarkers). Related disease folders are organized into a hierarchical tree. Gene content may vary greatly between such complex diseases as cancers and some Mendelian diseases. Also, coverage of different diseases in literature is skewed. These two factors may affect p-value prioritization for diseases; (d) Figure  legend giving the meanings of the colors and statistical interpretation of the above three figures. These figures were generated by MetaCore™. All results are sorted according to those most strongly represented in the "common" set. (a) Pathway maps. Canonical pathway maps represent a set of signaling and metabolic maps that comprehensively cover human biological functioning. All maps are created by Thomson Reuters scientists using a high quality manual curation process based on published peer reviewed literature. Experimental data are visualized on the maps as blue (for downregulation) and red (upregulation) histograms. The height of the histogram corresponds to the relative expression value for a particular gene/protein; (b) Process networks. The content of these cellular and molecular processes is defined and annotated by Thomson Reuters scientists. Each process represents a pre-set network of protein interactions characteristic of that process; (c) Related diseases (determined based on biomarkers). Related disease folders are organized into a hierarchical tree. Gene content may vary greatly between such complex diseases as cancers and some Mendelian diseases. Also, coverage of different diseases in literature is skewed. These two factors may affect p-value prioritization for diseases; (d) Figure legend giving the meanings of the colors and statistical interpretation of the above three figures. These figures were generated by MetaCore™.

The Highest-Scoring Pathway Map of Post-TBI-Related Biomarker Genes at the Four Time Points
The highest-scoring pathway map we identified using MetaCore for the four time points was DNA damage ATM (Ataxia-telangiectasia mutated )/ATR (Ataxia-telangiectasia and Rad3-related) regulation of G1/S checkpoint (Figure 3a). The ATM-Chk2 (ataxia telangiectasia mutated-checkpoint kinase 2) and ATR-Chk1 pathways are two distinct kinase signaling cascades that primarily coordinate cellular responses to DNA damage. The DNA damage checkpoints may arrest or delay progression of cell cycles in response to DNA damage. Notably, NF-κB, c-MYC, and P21 were up-regulated at all four time points relative to the time at which the injury occurred. Ubiquitin was down-regulated at both 4 and 8 h and up-regulated at 24 and 72 h. BRCA1 was down-regulated and MDM2 and I-κB were up-regulated at 4 h. PCNA was up-regulated at 24 h and CDK2 was up-regulated at 24 and 72 h. The second highest-scoring map was Cell Cycle ESR1 regulation of G1/S transition, which is also a cell cycle associated pathway (Figure 3b

The Highest-Scoring Pathway Map of Post-TBI-Related Biomarker Genes at the Four Time Points
The highest-scoring pathway map we identified using MetaCore for the four time points was DNA damage ATM (Ataxia-telangiectasia mutated )/ATR (Ataxia-telangiectasia and Rad3-related) regulation of G1/S checkpoint (Figure 3a). The ATM-Chk2 (ataxia telangiectasia mutated-checkpoint kinase 2) and ATR-Chk1 pathways are two distinct kinase signaling cascades that primarily coordinate cellular responses to DNA damage. The DNA damage checkpoints may arrest or delay progression of cell cycles in response to DNA damage. Notably, NF-κB, c-MYC, and P21 were up-regulated at all four time points relative to the time at which the injury occurred.     (Figure 4c) provide clues about new uses for conventional drugs and facilitate comparisons of the hidden molecular mechanisms common to some of these diseases.

The Highest-Scoring Pathway Map of Post-TBI-Related Biomarker Genes at 24 and 72 h
The highest-scoring pathway map at 24 and 72 h, according to the MetaCore™ analysis, is Cell Cycle Role of SCF Complex (Skp, Cullin, F-box containing complex) in Cell Cycle Regulation (Figure 5a). The S-phase kinase-associated protein (SKP1)/Cullin/F-box (SCF) complex is one of the E3-ubiquitin ligases, and plays a critical role in the cell cycle. Surprisingly, ubiquitin was down-regulated at 4 and 8 h, but up-regulated at 24 and 72 h. CDK1 and CDK2 were up-regulated at 24 and 72 h, and P21 was up-regulated at all four time points. Another pathway related to the cell cycle, influence of Ras and Rho proteins on G1/S transition, was also identified as strongly associated with biomarker genes at 24 and 72 h (Figure 5b). These genes, such as P21, CDK2, ERK1/2, MDM2, STAT3, RhoA, NF-κB, might provide clues for deciphering the regulatory mechanisms of post-TBI processes at different time points.    Results are sorted according to those most strongly represented in the "common" set. (a) Pathway maps. Canonical pathway maps represent a set of signaling and metabolic maps that comprehensively cover human biological functioning. All maps are created by Thomson Reuters scientists using a high quality manual curation process based on published peer reviewed literature. Experimental data is visualized on the maps as blue (for down-regulation) and red (up-regulation) histograms. The height of the histogram corresponds to the relative expression value for a particular gene/protein; (b) Gene ontology (GO) cellular processes. Since most GO processes have no gene/protein content, the "empty terms" are excluded from p-value calculations; (c) Process networks. The content of these cellular and molecular processes is defined and annotated by Thomson Reuters scientists. Each process represents a pre-set network of protein interactions characteristic for the process; (d) Figure legend giving the meanings of the colors and statistical interpretation of the above three figures. These figures were generated by MetaCore™. Canonical pathway maps represent a set of signaling and metabolic maps that comprehensively cover human biological functioning. All maps are created by Thomson Reuters scientists using a high quality manual curation process based on published peer reviewed literature. Experimental data is visualized on the maps as blue (for down-regulation) and red (up-regulation) histograms. The height of the histogram corresponds to the relative expression value for a particular gene/protein; (b) Gene ontology (GO) cellular processes. Since most GO processes have no gene/protein content, the "empty terms" are excluded from p-value calculations; (c) Process networks. The content of these cellular and molecular processes is defined and annotated by Thomson Reuters scientists. Each process represents a pre-set network of protein interactions characteristic for the process; (d) Figure legend

Top-Scoring AN (Analyze Networks Algorithm) Results for the Four Time Points Post-TBI
The AN results for the four time points provide information about canonical pathways, up-regulated and down-regulated genes, and mixed expression genes ( Figure 6).

Discussion of the Cell Cycle Behavior of TBI
Our results show that the cell cycle is strongly related to brain injury (Figure 1). Wu et al. [30] discuss cell cycle activation and cord injury, stating that a traumatic spinal cord injury (SCI) causes a series of events involving initial mechanical damage, secondary injury processes, and eventually results in tissue loss and functional impairment. They also found that the cell cycle is activated following an SCI [30].

Network Ontology Analysis (NOA) Results
The detailed analytical results of the network ontology analysis (NOA) highlight the evolutionary process at the four time points post-TBI with respect to biological processes, cellular components, and molecular functions (Table 4). For example, the most statistically significant biological process at 4 h post-TBI was protein modification by small-protein conjugation, while the top processes at the other three time points were positive regulation of cellular metabolic processes. The top cellular components at the four time points were protein complex, nucleus, intracellular part, and cytoplasm, respectively. The top molecular function at all four time points was protein binding, followed by enzyme binding, promoter binding, receptor signaling protein activity, and identical protein binding, at 4, 8, 24, and 72 h post-TBI, respectively. The AN results for the four time points provide information about canonical pathways, up-regulated and down-regulated genes, and mixed expression genes ( Figure 6).

Discussion of the Cell Cycle Behavior of TBI
Our results show that the cell cycle is strongly related to brain injury (Figure 1). Wu et al. [30] discuss cell cycle activation and cord injury, stating that a traumatic spinal cord injury (SCI) causes a series of events involving initial mechanical damage, secondary injury processes, and eventually results in tissue loss and functional impairment. They also found that the cell cycle is activated following an SCI [30].

Network Ontology Analysis (NOA) Results
The detailed analytical results of the network ontology analysis (NOA) highlight the evolutionary process at the four time points post-TBI with respect to biological processes, cellular components, and molecular functions (Table 4). For example, the most statistically significant biological process at 4 h post-TBI was protein modification by small-protein conjugation, while the top processes at the other three time points were positive regulation of cellular metabolic processes. The top cellular components at the four time points were protein complex, nucleus, intracellular part, and cytoplasm, respectively. The top molecular function at all four time points was protein binding, followed by enzyme binding, promoter binding, receptor signaling protein activity, and identical protein binding, at 4, 8, 24, and 72 h post-TBI, respectively.
(a)   Up-regulated genes are marked with red circles, and down-regulated with blue circles. A "checkerboard" pattern indicates mixed expression for the gene between files or between multiple tags for the same gene. These figures were generated by MetaCore™. Table 4 In addition to our original analysis that identified the PPINs and their corresponding network biomarkers, we have generated abundant results using DAVID (The Database for Annotation, Visualization and Integrated Discovery; Available online: https://david.ncifcrf.gov/home.jsp), MetaCore™, and NOA. Based on our previous experience with cancer research, we applied the David analysis to the four time points post-TBI and expected it to yield the results of KEGG (Kyoto Encyclopedia of Genes and Genomes; Available online: http://www.genome.jp/kegg/). However, we did not find the KEGG database suitable for analyzing TBI (most entries were linked to cancer-related pathways), so we present these results in Supplementary Material S4 and do not describe them here.

Summary of the Results in
Here, we emphasize the importance of the results of the MetaCore™ and NOA analyses. They can be used as references for medical and basic researchers performing further explorations of the hidden mechanisms of network biomarkers, and will also help these researchers to develop novel strategies for TBI therapy or recovery processes. Our original results can also be used for further analysis using other methods, so we present them here.  Table 4 In addition to our original analysis that identified the PPINs and their corresponding network biomarkers, we have generated abundant results using DAVID (The Database for Annotation, Visualization and Integrated Discovery; Available online: https://david.ncifcrf.gov/home.jsp), MetaCore™, and NOA. Based on our previous experience with cancer research, we applied the David analysis to the four time points post-TBI and expected it to yield the results of KEGG (Kyoto Encyclopedia of Genes and Genomes; Available online: http://www.genome.jp/kegg/). However, we did not find the KEGG database suitable for analyzing TBI (most entries were linked to cancer-related pathways), so we present these results in Supplementary Material S4 and do not describe them here.

Summary of the Results in
Here, we emphasize the importance of the results of the MetaCore™ and NOA analyses. They can be used as references for medical and basic researchers performing further explorations of the hidden mechanisms of network biomarkers, and will also help these researchers to develop novel strategies for TBI therapy or recovery processes. Our original results can also be used for further analysis using other methods, so we present them here. Table 4. Network ontology pathway analysis and gene set enrichment analysis for four time points post-TBI with respect to (1) biological processes; (2) cellular components; and (3) molecular functions. Each of these tables is presented for each of the four time points. R: number of genes in reference set; T: number of genes in test set; G: number of genes annotated by given term in reference set; O: number of genes annotated by given term in test set. This table is produced based on data stored in the network ontology analysis (NOA) web server.

Overview of the Traumatic Brain Injury (TBI) Network Biomarkers Identification Process
We used a theoretical framework to identify the evolution of TBI network biomarkers at four time points representing four important stages after the occurrence of a TBI. We have successfully used a similar theoretical framework to identify the network biomarkers of various cancers and stroke [31][32][33]. We therefore highlight only the key points of this framework in this section, and describe the process in detail in the supplementary information (Supplementary Material S0). The flowchart in Figure 7 illustrates how we identified the network biomarkers of TBI at four time points. We combined two kinds of data sets to develop our model: (i) microarray data from both TBI and normal samples from the Gene Expression Omnibus (GEO) database, where TBI samples are divided into four groups according to the time post-injury (4,8,24, and 72 h); and (ii) a protein-protein interaction database, which was necessary to construct four candidate PPINs for TBI.

Overview of the Traumatic Brain Injury (TBI) Network Biomarkers Identification Process
We used a theoretical framework to identify the evolution of TBI network biomarkers at four time points representing four important stages after the occurrence of a TBI. We have successfully used a similar theoretical framework to identify the network biomarkers of various cancers and stroke [31][32][33]. We therefore highlight only the key points of this framework in this section, and describe the process in detail in the supplementary information (Supplementary Material S0). The flowchart in Figure 7 illustrates how we identified the network biomarkers of TBI at four time points. We combined two kinds of data sets to develop our model: (i) microarray data from both TBI and normal samples from the Gene Expression Omnibus (GEO) database, where TBI samples are divided into four groups according to the time post-injury (4,8,24, and 72 h); and (ii) a protein-protein interaction database, which was necessary to construct four candidate PPINs for TBI.  Using several mathematical models, we constructed four TBI PPINs (TPPINs) and a normal PPIN (NPPIN). We also constructed four Differential PPINs (DPPINs) and calculated the TBI relevance value (TRV) for each protein in the DPPINs (Supplementary Material S0).

Data Selection and Pre-Processing
To identify the regulatory mechanisms of brain trauma, we obtained microarray datasets from gene expression profiling of brain injury samples (4,8,24, and 72 h after surgery) and brain sham samples (4,8,24, and 72 h after sham surgery) in mice from the NCBI GEO [34]. We chose GSE2392 [35] and its corresponding platform, GPL81, as our research object (Table 1). We collected expression profiling data for mice from the Gene Expression Omnibus (GEO), accession number GSE2392, and normalized it using quantile normalization. For each of the four time points, there were three TBI samples and two control samples, making a total of 20 samples.
The PPI dataset for Mice was downloaded from the public Biological General Repository for Interaction Database (BioGRID) [36]. We excluded false-positive PPIs from the TBI PPINs using BioGRID and compared the resulting PPINs to identify network biomarkers (Supplementary Material S0).

Selection of the Differential Protein Groups and Identification of PPINs
We combined the gene expression profiles with the PPI database to construct the corresponding TPPINs and NPPIN. First, we identified a differential protein group containing proteins with differential expression. Protein expression profiles are not yet readily available, so we assumed that gene expression profiles are correlated with protein expression profiles. We used a one-way analysis of variance (ANOVA) or the fold change (FC) method to analyze the expression of each protein and select the proteins that were differentially expressed in cells. These methods were capable of identifying the critical proteins differentially expressed in TBI and normal samples. We eliminated proteins with no PPI information from the differential protein group. Proteins that were not already in the differential protein group were included if their PPI information indicated that they had a close relationship with a protein (number of edges > 5) in the group. Thus, the differential protein group included the critical differential expression proteins and other closely related proteins.
We combined the critical differential protein group and PPI information to construct candidate PPINs by linking proteins that interacted with each other. In other words, proteins for which we had PPI information via the differential protein group were linked using our mathematical model to form the candidate TPPINs and NPPINs.
The candidate TPPINs and NPPINs included all possible PPIs under various biological conditions. It was then necessary to confirm them using microarray expression data and identify the proper PPIs according to the TBI processes. We used both the PPI model selection method and a model order detection strategy to exclude the false-positive PPIs. We could then express the relationship to the PPI of the i-th target protein in the candidate TPPINs and NPPINs with the following formula: where x i (n) is the expression levels of target protein i for sample n; x j (n) is the expression level of the j-th protein interacting with target protein i for sample n; α ij is the association ability between target protein i and its j-th interactive protein; M i is the number of proteins interacting with target protein i; and ω i (n) is stochastic noise due to other factors or model uncertainty. For the definitions of the terms, please refer to Supplementary Material S1.
We then used the least squares parameter estimation method [37] to determine the parameter of association strength in Equation (1) from the TBI microarray data as follows (Supplementary Material S1): x i pnq " M i ÿ j"1α ij x j pnq, i " 1, 2,¨¨¨,M whereα ij is solved using TBI microarray data and the least squares parameter estimation method. Finally, we use AIC (Akaike information criterion) [37] with a Student's t-test [38] for model order determination and to determine the critical protein association strengthsα ij (Supplementary Material S2).

Determination of Proteins with Top TBI Relevance Value (TRV)and Their Corresponding Network Structures
The significant remaining PPIs after pruning the false-positives are expressed by the following: where M i 1 ď M i is the total number of significant PPIs remaining in the refined PPINs for target protein i. The final refined PPIN is as follows: Xpnq " AXpnq`wpnq (4) Xpnq " » ----- . .  The interaction matrix A of refined PPINs (i.e., after pruning, as described by Equation (4) were constructed as follows: where k is time post-injury, A K T and A N are the interaction strength matrices of the refined TPPINs and NPPIN for 4, 8, 24, and 72 h post-injury, respectively, and M is the total number of interaction proteins in the refined TPPINs and NPPIN after pruning. Thus, the protein interaction model can be expressed as follows: where k is time post-injury and x k T pnq " rx k 1T x k 2T¨¨¨x k MT s T and x N pnq " rx 1N x 2N¨¨¨xMN s T are vectors of protein expression values. The difference matrix D k we used to interpret the differences in behavior between two networks is defined as follows: where k is time post-injury, d k ij is the difference between the protein association capacity of the networks at different times post-injury, and matrix D k denotes the differences in network structures. To investigate TBI-related factors using matrix D k , we proposed a value we called the TBI relevance value (TRV) to quantify the significance of each protein in D k as follows [13]:

Pathway Analysis
Our theory allowed us to construct the TPPINs and NPPINs and identify the network biomarkers with the highest TRVs. However, we also wanted to explore the biological significance of these network biomarkers; here, we used various software packages and on-line tools to do the functional pathway analysis. For example, we used KEGG [39], the DAVID bioinformatics database [40,41], and NOA [42,43], which are the most commonly used free on-line tools, and MetaCore™ from Thomson Reuters, a powerful commercial software package that can do multiple functional and pathway analyses. A detailed description of the process we followed is given in Supplementary Material S0.

Software and Databases Used in This Research
analysis; Available online: http://apps.cytoscape.org/apps/noa) are easy-to-use web servers with user-friendly interfaces.

Conclusions
In conclusion, we used a publicly available high-throughput microarray to construct networks for four time points post-TBI. The powerful functional pathway analysis then extended the relevance of our research to a wider field. This work can help scientists and medical researchers reveal the underlying hidden mechanisms of TBI and improve biological models for TBI. However, the number of samples in this study was small. In future, we will include more data samples from humans and model organisms to develop a more precise and accurate model. In addition, we will integrate the GRN and PPI models to facilitate deeper interpretation of the regulatory and interaction mechanisms. This work can offer clues and information for developing novel biological models and strategies for target therapy to advance recovery processes.