Graph theory analysis of DTI tractography in children with traumatic injury.

OBJECTIVE: To evaluate brai structural connectivity in children with traumatic injury (TI) following a motor vehicle accident using graph theory analysis of DTI tractography data.
METHODS: DTI scans were acquired on a 3 T Philips scanner from children aged 8-15 years approximately 2 months post-injury. The TI group consisted of children with traumatic brain injury (TBI; n = 44) or extracranial injury (EI; n = 23). Healthy control children (n = 36) were included as an age-matched comparison group. A graph theory approach was applied to DTI tractography data to investigate injury-related differences in connectivity network characteristics. Group differences in structural connectivity evidenced by graph metrics including efficiency, strength, and modularity were assessed using the multi-threshold permutation correction (MTPC) and network-based statistic (NBS) methods.
RESULTS: At the global network level, global efficiency and mean network strength were lower, and modularity was higher, in the TBI than in the control group. Similarly, strength was lower and modularity higher when comparing the EI to the control group. At the vertex level, nodal efficiency, vertex strength, and average shortest path length were different between all pairwise comparisons of the three groups. Both nodal efficiency and vertex strength were higher in the control than in the EI group, which in turn were higher than in the TBI group. The opposite between-group relationships were seen with path length. These between-group differences were distributed throughout the brain, in both hemispheres. NBS analysis resulted in a cluster of 22 regions and 21 edges with significantly lower connectivity in the TBI group compared to controls. This cluster predominantly involves the frontal lobe and subcortical gray matter structures in both hemispheres.
CONCLUSIONS: Graph theory analysis of DTI tractography showed diffuse differences in structural brain network connectivity in children 2 months post-TI. Network differences were consistent with lower network integration and higher segregation in the injured groups compared to healthy controls. Findings suggest that inclusion of trauma-exposed comparison groups in studies of TBI outcome is warranted to better characterize the indirect effect of stress on brain networks.

Introduction

Pediatric traumatic injury (TI) is a major public health problem. >9 million children—about 1 in 4—visit the emergency department for treatment of injuries each year in the United States (Borse and Sleet, 2009). Traumatic brain injury (TBI), a subset of TI, is the most common cause of death and acquired brain injury among youth in the United States. TBI occurs in >100,000 children each year and incurs an estimated annual cost of more than $1 billion (Graves et al., 2015). After TI due to motor vehicle collision in particular, children with TBI and extracranial bodily injury (EI) are at risk of experiencing traumatic stress (Max et al., 1998). Indeed, 1 in 3 survivors of TI are estimated to develop posttraumatic stress disorder, and even more will experience significant posttraumatic stress symptoms including hyperarousal, avoidance, re-experiencing, and emotional numbing (Holbrook et al., 2005; Langeland and Olff, 2008). Other commonly occurring co-morbid conditions following TI include anxiety, depression, and substance abuse, as well as reduced adaptive behavior and health-related quality of life (Fay et al., 2009; Max et al., 2012; Zatzick et al., 2008). Such psychological health problems, of which TI is a significant risk factor, can far outlast the physical injury.

Although there is a substantial literature regarding changes to brain structure following TBI, very little is known regarding the impact of bodily injury on brain structure in other forms of TI. In TBI, multifocal tissue damage occurs at macrostructural and microstructural levels resulting from rotational and translational forces. Frontal and temporal regions, and especially the hippocampus, are most susceptible to injury (DeMaster et al., 2017; Wilde et al., 2012; Wilde et al., 2007; Wilde et al., 2005). Beyond direct injury to the brain following impact, it is possible that the psychological health problems after TI might also result from the stress response disrupting functioning of prefrontal cortex and the limbic system, particularly the amygdala and hippocampus (Carrion et al., 2007; Jackowski et al., 2009; Juranek et al., 2012; Karl et al., 2006). These brain regions continue to mature through adolescence and, being central to the hypothalamic-pituitary-adrenal (HPA) axis, are vulnerable to increased stress (McEwen et al., 2015; McEwen et al., 2016). Consequently, limbic-prefrontal structures (i.e., prefrontal cortex, amygdala, hippocampus, and anterior cingulate cortex) may be disproportionately affected by injury to the brain or other body regions, as well as by traumatic stress. Furthermore, these regions are especially likely to show aberrant connectivity when extreme stress or mechanical brain injury occur in childhood or adolescence (Nooner et al., 2013).

Diffusion tensor imaging (DTI) has been employed in many TBI studies due to its sensitivity to white matter (WM) microstructure. During the past decade, a number of DTI studies of pediatric TBI have revealed widespread injury to commissural, association, and projection pathways (Ewing-Cobbs et al., 2016; Genc et al., 2017; Ryan et al., 2018). While WM microstructure has been extensively studied in pediatric TBI, very little research using DTI has been carried out in non-specific TI or in EI. There is, however, some research in children who have experienced significant stress in the context of posttraumatic stress disorder (PTSD), often related to maltreatment. The majority of analyses have been restricted to a few major WM tracts, but tend to show reductions in fractional anisotropy (FA) in the corpus callosum, uncinate fasciculus, and inferior and superior longitudinal fasciculi after childhood trauma (Choi et al., 2012; Choi et al., 2009; Jackowski et al., 2008). In the months following injury, studies consistently find that FA is lower and MD is higher in children with TBI, particularly those with moderate to severe injury (Ewing-Cobbs et al., 2016; Genc et al., 2017). These differences are considered to reflect WM injury and the impact of ongoing neurodegenerative processes (DeKosky and Asken, 2017). Importantly, this WM injury is present in the majority of tracts that have been studied, highlighting the global effects of TBI on brain connectivity. To our knowledge, there have been no DTI studies of children with TI as the primary focus, except for inclusion of children with non-brain injuries as “trauma controls” to compare to TBI. There is limited understanding of how TI alters the structural connectivity between distributed brain networks. Moreover, it is unclear how injury sustained during childhood or adolescence may affect the structural architecture of brain circuitry with protracted developmental trajectories such as the prefrontal-limbic and fronto-parietal networks. Due to the substantial development of WM connectivity during childhood and adolescence, TI experienced during these developmental stages may significantly disrupt functioning during the early stages of recovery and may have a disproportionately negative impact on subsequent brain maturation and strengthening of connectivity of key neural circuitry.

Application of graph analysis to DTI data may illuminate the diffuse consequences of TI by characterizing changes in brain network organization (Irimia et al., 2012). In the graph-theoretical approach, the brain is represented as a network (or graph in mathematics) consisting of vertices (brain regions or nodes) and edges (connections between regions) (Newman, 2010). With this representation, a number of metrics can be calculated which are measures of different aspects of global or local network connectivity. Characteristic path length is a measure of integration of information processing; an increase in this metric implies that it takes more “steps” for any two brain regions, on average, to communicate with one another. Global efficiency, another measure of integration, quantifies a network's parallel information processing efficiency in a network and is proportional to the inverse path length (Latora and Marchiori, 2001). Strength is simply a weighted average of the measured connectivity strength among all of a region's connections. Finally, modularity is a measure of network segregation, calculated by partitioning a network into groups of regions (modules, or communities) with high connectivity within modules relative to the connectivity between regions in distinct modules (Girvan and Newman, 2002; Newman and Girvan, 2004). Segregation is related to the concept of functional specialization (Baum et al., 2017). Metrics from healthy adult brains typically have lower characteristic path length, higher global efficiency, higher strength, and lower modularity. Across childhood and adolescence, developmental changes reflecting greater network integration are associated with increases in global efficiency and vertex/node strength and decreases in path length and modularity (Dennis et al., 2013; Hagmann et al., 2010).

To date, only a few studies have used graph analysis to investigate brain structural connectivity in children after TBI, with discrepant findings which may be due to differences in: patient groups (i.e., injury severity); comparison groups (healthy or TI); time since injury; and network construction methods. Caeyenberghs et al., 2012 created networks from DTI tractography using “fiber count” as a measure of connectivity strength in adolescents with moderate-to-severe TBI approximately 3 years after injury. They reported that characteristic path length was increased and local efficiency was decreased across the brain, compared to a healthy control group (Caeyenberghs et al., 2012). More recently, increased modularity and path length was also evident using fiber count-based tractography networks for adolescents with acute mild TBI (mTBI) compared to children with orthopedic injury (Yuan et al., 2015). In contrast to Yuan and colleagues, Konigs et al., 2017 did not find differences in any network-based metrics between children with mTBI and TI. In agreement with Caeyenberghs et al., 2012 who also acquired DTI 3 years post-injury, characteristic path length was increased in the moderate-to-severe TBI group (Konigs et al., 2017). More recently, Dennis et al. studied children and adolescents with post-acute moderate-to-severe TBI, and found network differences between children who experienced post-traumatic seizures and those who did not (Dennis et al., 2017). Although the limited available studies provide important information regarding the effects of TBI on brain network structure in children, little is known about the direct or indirect impact of extracranial injuries. Moreover, there is a dearth of information regarding the effects of different levels of injury severity or the effects of demographic variables such as age, sex, or pubertal status on brain architecture after injury.

To address the gaps in the literature, our aim is to examine the impact of microstructural changes post-TI on brain connectivity. We used a graph theoretical approach to analyze brain WM networks in children approximately 2 months after injury to determine how changes in the integrity of the brain structural network differ from those in children with injuries to other body regions as well as in healthy controls. To do so, we used the multi-threshold permutation correction (MTPC) approach (Drakesmith et al., 2015), which assesses differences in either global or local network metrics, and the network-based statistic (NBS) (Zalesky et al., 2010), which assesses differences in the strength of each individual connection. Our central hypothesis is that TI disrupts—directly and/or indirectly—brain network connectivity globally and between prefrontal and limbic structures specifically. We expected both EI and TBI sub-groups to differ from the control group, with the TBI group showing greater deviation in network metrics from the control group than the EI group. In addition, network metrics were expected to differ based on severity of TI, such that more severe injury is associated with more disrupted networks, and age, such that older age at injury is associated with less disrupted networks.

Methods

Subjects

This study is part of a prospective longitudinal study. Youth ages 8–15 years sustaining TBI or EI in motor vehicular-related accidents (as passengers or pedestrians) were recruited from the Emergency Department or Level 1 Pediatric Trauma Center at Children's Memorial Hermann Hospital/University of Texas Health Science Center at Houston (UTHealth) between September 2011 and August 2016. Participants with TBI and EI met the following inclusion criteria: 1) injured in a motor vehicle accident between ages 8 and 15 years (further described in the next paragraph); 2) proficiency in English or Spanish; 3) residing within a 125 mile catchment radius; 4) no prior history of major neuropsychiatric disorder (intellectual deficiency or low-functioning autism spectrum disorder) that would complicate assessment of the impact of injury on brain outcomes; 5) no metabolic, endocrine, or systemic health problems (e.g., hypertension); 6) no prior medically-attended TBI; and 7) no habitual use of steroids, tobacco, or alcohol. The latter four criteria were assessed during screening using a brief parent interview.

For the present study, “motor vehicular-related accidents as passengers or pedestrians” includes the involvement of any type of motorized vehicle, defined as a vehicle not operating on rail and powered by a motor. Regarding the “mechanism of injury” in Table 1, “bicycle” injuries were those in which the subject was riding a bicycle and struck by a motor vehicle. Furthermore, “fall from vehicle” injuries occurred when the vehicle, such as a Quad 4Wheeler, was in motion.

Table 1

Subject demographic and injury variables.

Values are N (%) or median (IQR). P-values were calculated using the Kruskal-Wallis rank sum test for continuous variables, and Fisher's exact test for categorical variables. P-values for injury-related variables were calculated based on comparison between the EI and TBI groups only, as the control group did not have values for those variables.

	Control (N = 36)	Extracranial Injury (N = 21)	Traumatic Brain Injury (N = 44)	P-value
Demographic Variables
Age at MRI, mo	144 (128–174)	150 (119–172)	153 (131–172)	0.85
Sex, F	15 (42)	8 (38)	19 (43)	0.97
Scanner change, Post	13 (36)	5 (24)	18 (41)	0.42
Pubertal Development Scale	3 (2–3)	2 (1–3)	3 (2–3)	0.93
Psychosocial Adversity Index	1 (1–2)	1 (1–2)	1 (1–2)	0.97
Ethnicity, Hispanic	24 (67)	12 (57)	23 (52)	0.44
Race				0.69
White	25 (69)	17 (81)	30 (68)	–
African American	8 (22)	4 (19)	12 (27)	–
Asian	0 (0)	0 (0)	1 (2)	–
Multi	3 (8)	0 (0)	1 (2)	–
Injury Variables
Mechanism of injury				0.39
Bicycle	–	8 (38)	16 (36)	–
Motor vehicle collision	–	4 (19)	14 (32)	–
All-terrain vehicle	–	4 (19)	7 (16)	–
Motorcycle	–	3 (14)	1 (2)	–
Fall from vehicle	–	2 (10)	6 (14)	–
Time since injury, wk	–	7 (6–9)	7 (6–8)	0.17
Injury Severity Score (ISS)	–	10 (5–14)	11 (9–17)	0.24
ISS (no head)	–	10 (5–14)	2 (1–5)	<0.001
Abbreviated Injury Scale ≥ 2
Head and neck	–	1 (5)	38 (86)	<0.001
Face	–	0 (0)	3 (7)	0.55
Chest	–	6 (29)	11 (25)	0.77
Abdomen	–	4 (19)	0 (0)	0.009
Extremity	–	18 (86)	9 (20)	<0.001
External	–	0 (0)	4 (9)	0.30
Lowest Glasgow Coma Scale
3–8 (severe)	–	–	11 (25)	–
9–12 (moderate)	–	–	4 (9)	–
13–15 (mild)	–	–	29 (66)	–

mo: months; F: female; wk: weeks; ISS: Injury Severity Score; GCS: Glasgow Coma Scale.

Of the 220 injured youth who met study inclusion criteria, 131 were consented and enrolled, and 112 were scanned at the 2 month follow-up (TBI: n = 80; EI: n = 32). Youth with EI were included to investigate the effects of different types of traumatic injuries on brain network connectivity while also allowing us to account for stresses of injury that may influence brain development. We removed 2 subjects from the EI group due to their having history of traumatic (extracranial) injury. A healthy comparison group was recruited from the community and met criteria 2–7; of 78 eligible youth, 56 were consented and 52 were scanned. Informed written consent was obtained from the guardian of each participant and written assent was obtained from all participants in accordance with Institutional Review Board guidelines.

Procedure

Exposure to adversity prior to the injury was assessed by trained interviewers at the time of enrollment and included the following common indicators: 1) severe marital discord defined as parental divorce or separation; 2) low social status defined as levels IV or V on the Hollingshead Index; 3) large family size defined as three or more children living in the child's primary home; 4) history of investigation by protective service agency regarding this child; 5) parental criminal conviction; and 6) treatment of parental mental health problems (Biederman et al., 2002). Items were scored “yes/no” and summed to yield a total Psychosocial Adversity Index (PAI).

Pubertal changes, including growth, body hair, skin, and sex-specific changes were independently rated by both children and their parents at the time of scanning using the Petersen Pubertal Development Scale (PDS) (Petersen et al., 1988). Each item was then coded on a 5-point scale similar to Tanner staging (Shirtcliff and Essex, 2008); ratings were averaged to yield a score ranging from 1 (pre-pubertal) to 5 (post-pubertal). If ratings of parents and children differed by > 1 point, they were asked to discuss and come to consensus, and the consensus rating was used in analyses.

Subjects in the TBI group were further subdivided based on injury severity using the Glasgow Coma Scale (GCS) score, obtained upon admission to the hospital (Teasdale and Jennett, 1974). The GCS consists of a 1–5 rating in three domains: eye responses (eye opening to stimuli), verbal response, and motor response. The total score is a sum of the three domains and ranges from 3 (coma or death) to 15 (normal). We defined 3 subgroups: mild (GCS 13–15), complicated-mild (GCS 13–15 combined with acute hemorrhage or parenchymal injury), and moderate-severe (GCS 3–12). Injury severity in the EI and TBI groups was measured using the Abbreviated Injury Scale (AIS), which classifies severity of injury to specified anatomical regions (head/neck, face, chest, abdomen, extremities, and external/other) on a scale from 1 to 6. The injury severity score (ISS) is the sum of the square of the highest AIS scores from three anatomical regions and ranges from 0 to 75 (S. P. Baker et al., 1974). We calculated both the total ISS and an ISS excluding head/neck injury. We report the number of subjects with AIS ≥ 2, excluding subjects with minor injury, for each anatomical region in Table 1.

MRI acquisition

Subjects were scanned on a Philips 3 Tesla (T) Intera system (Philips Healthcare, Amsterdam, Netherlands) with a 32-channel head coil at University of Texas McGovern Medical School. During the study, several datasets were corrupted due to scanner equipment failure (n = 8); these subjects were excluded from the present study. In November 2014, the scanner was upgraded to a Philips 3 T Ingenia; as such, we adjust for scanner hardware upgrade in our statistical analyses. The T1-weighted sequence was acquired in the sagittal plane with parameters: TR/TE = 8.07/3.68 ms, flip angle = 6°, acquisition matrix = 256 × 256, FOV = 256 mm, slice thickness = 1 mm, with resultant voxel size = 1 × 1 × 1 mm3. The DTI sequence was acquired using single-shot spin-echo echo planar imaging with parameters: TR/TE = 8700/67 ms, flip angle = 90°, acquisition matrix = 96 × 96, FOV = 240 mm, slice thickness = 2.5 mm, with resultant voxel size = 2.5 × 2.5 × 2.5 mm3. A single non-diffusion weighted volume was acquired (b = 0 s / mm2), along with 32 diffusion-weighted volumes (b = 1000 s / mm2).

DTI preprocessing and tractography

To define regions of interest (ROI), each subject's T1-weighted image was parcellated using Freesurfer v5.3.0 and the Desikan-Killiany (DK) atlas (Dale et al., 1999; Desikan et al., 2006; Fischl and Dale, 2000; Fischl et al., 1999). Surfaces and parcellations were manually inspected and corrected if errors were present. The 82 regions of the DK atlas (68 cortical and 14 subcortical) were registered to diffusion space using boundary-based registration (BBR) in Freesurfer (Greve and Fischl, 2009). In addition to the 82 gray matter regions, a mask of the ventricles was created and registered to diffusion space. All registrations were checked visually for accuracy. Subjects with large lesions which could not be processed accurately (N = 1) were excluded from further analyses. Image quality was assessed visually by 2 of the study's authors (CGW and DD), and subjects were excluded based on agreement between the two.

All DTI image processing was performed with FSL v5.0.9 and the FSL Diffusion Toolbox (FDT) (Smith et al., 2004). We used the Texas Advanced Computing Center (TACC) (University of Texas) for all DTI processing steps. First, eddy-current correction was performed using the latest GPU version of eddy to correct for image distortions and head motion (Andersson et al., 2016; Andersson et al., 2003; Andersson and Sotiropoulos, 2015, Andersson and Sotiropoulos, 2016; Graham et al., 2016). Distributions of diffusion parameters were calculated at every voxel using the latest GPU version of bedpostx with default parameters (Behrens et al., 2007; Behrens et al., 2003; Hernandez et al., 2013; Jbabdi et al., 2007). This program fits parameters for 3 fibers per voxel with the “ball and multi-stick with a Gamma distribution of diffusivities” deconvolution model (Jbabdi et al., 2012).

Probabilistic tractography was performed with the tool probtrackx2 (Behrens et al., 2007, 2003). Fiber tracking was initiated from each seed region individually, with the remaining 81 regions as targets. For each seed voxel, 5000 samples were drawn from the probability distribution of the principal fiber direction (estimated by bedpostx); each step in the process was 0.5 mm, and tracking stopped after a maximum of 2000 steps. The ventricles mask was specified as an “avoid” mask, so that streamlines entering the ventricles were rejected. Additionally, we corrected the resulting path distributions for pathway lengths: for each seed region, the average length of all samples reaching a target was calculated and multiplied by the number of samples reaching that target. All other parameters used were program defaults.

Network construction

Network construction and statistical analyses were performed in R version 3.4.3 (2017-11-30), using the packages igraph v1.2.1 and brainGraph v2.4.0 (Csardi and Nepusz, 2006; Kolaczyk and Csardi, 2014; R Core Team, 2015; Watson, 2018). For each subject, a connectivity matrix A was obtained from the probabilistic tractography procedure described above. The matrices for all subjects were combined into a multidimensional array, such that each element A(i, j, k) equals the number of streamlines (i.e., fibers that connect two regions), weighted by pathway length, between regions of interest (ROIs) i and j for subject k. The operations performed on these matrices to remove spurious connections are described in the remainder of this section.

Pairs of larger brain regions have higher streamline counts by virtue of being larger targets. To account for the variation in ROI size, we divided each matrix entry by the average volume of the ROI pairs (Gong et al., 2009; Hagmann et al., 2008). Next, a set of thresholds T (n = 30) was applied at the individual level to remove the lowest connectivity weights (i.e., entries less than that threshold were changed to 0). These values were chosen empirically and resulted in a range of densities similar to the literature (range: 4.9–31.1%) (Gong et al., 2009; Hagmann et al., 2008).

An additional threshold was applied at the individual subject level to balance the number of false positive and false negative structural connections. Conducted separately for each group (i.e., TBI, EI, and control), only connections that were present in at least 50% of subjects per group were retained. This group threshold (50%) is within the optimal range determined by a previous report (de Reus and van den Heuvel, 2013). In sum, 30 connectivity matrices per subject were generated through this thresholding process. The final edge weights are pathway length-corrected streamline counts normalized by average ROI-pair volumes.

FA-weighted networks

After we generated networks weighted by normalized streamline counts, we calculated mean FA along the streamlines between ROI pairs and set these values as edge weights of the networks (Konigs et al., 2017). To remove the influence of spurious streamlines when calculating mean FA, the path distribution was thresholded such that only the top 10% voxels (in terms of number of streamlines) were included when calculating mean FA.

Statistical analysis

Between-group differences in demographic variables were assessed using the Kruskal-Wallis rank sum test for continuous variables, and Fisher's exact test for categorical variables. Differences in injury-related variables were calculated between the EI and TBI groups only. For between-group comparisons of network metrics, we specified General Linear Models (GLMs) adjusting for age at MRI, sex, and scanner (either pre- or post-hardware update). From this model, we added interaction terms to test for significant interactions: group X sex, group X age, and group X scanner, but removed the interaction terms from the final models due to lack of statistical significance for all network measures tested. Analyses of the associations between network metrics and injury severity included both the EI and TBI groups in a single GLM; a group X ISS interaction was tested but removed from the final model due to lack of statistical significance.

Network metrics

Between-group differences in network metrics were assessed using the multi-threshold permutation correction (MTPC) method (Drakesmith et al., 2015). This approach determines significance based on the presence of significant effects across multiple adjacent thresholds. MTPC was applied to both vertex- (strength, nodal efficiency, local efficiency, average shortest path length, and weighted participation coefficient) and graph-level metrics (global efficiency, characteristic path length, modularity, and strength). Briefly, a GLM is specified at each threshold (and, for vertex-level metrics, at each vertex) to test for a significant group difference in the network metric. The data were permuted according to the Freedman-Lane procedure (Freedman and Lane, 1983; Winkler et al., 2014). The same GLM is tested on the permuted data (at each vertex, if applicable) and the maximum statistic is recorded, building a null distribution.

We repeated this procedure 10,000 times for graph-level and 5,000 times for vertex-level measures per threshold. A critical statistical threshold (Scrit) is calculated as the 1 − αth percentile of the set of all permuted (maximum) statistics (here, α = 0.05). The area under the curve (AUC) is calculated for the observed test statistics S where the observed statistic exceeded Scrit for 3 or more consecutive thresholds (denoted Amtpc),

Similarly, a critical AUC (Acrit) was calculated as the average of the AUC's for the null statistics (under the same conditions as the calculation of Amtpc). A significant difference is present if Amtpc > Acrit.

Network-based statistic

To determine group differences in connection-wise connectivity strength, the network-based statistic (NBS) method was employed (Zalesky et al., 2010). This method allows for family-wise error control of network data that is analogous to cluster-based thresholding in the functional MRI literature. To calculate the NBS, first a GLM adjusting for age at MRI, sex, and scanner (either pre- or post-hardware update) was specified for each element of the N × N connectivity matrix (i.e., the matrices thresholded by ROI size). A N × N matrix of t-statistics associated with the contrast of interest was thresholded by an initial P-value threshold (here, P < .001). A graph was then created from this matrix, and the largest connected component recorded. Next, the data were permuted 5,000 times in which each subject was randomly assigned to one of the subject groups (of equal size as the original control, EI, and TBI groups). The same GLM was again specified at every matrix entry, a t-statistic matrix calculated for the permuted dataset, and the associated P-values were thresholded (again by P < .001). Finally, the largest connected component was recorded for the resultant graph; this procedure was repeated for each permutation. The null distribution of largest connected component sizes was used to calculate a P-value associated with the connected components of the observed data (as the proportion of times the permutation component sizes exceeded the observed size).

Results

A flowchart of subject inclusion and exclusion numbers is shown in Supplementary Fig. 1 (S1). Of the 164 subjects scanned, 63 were excluded; the most common reason was motion-related artifact (the DTI sequence; N = 37). Proportionally, more TBI subjects were excluded (28%) compared to control (17%) or EI subjects (16%). Furthermore, a comparison of included vs. excluded subjects found a significant difference only in ethnicity, such that a higher percentage of subjects included were of Hispanic ethnicity (P = .025). A total of 44 TBI subjects, 21 EI subjects, and 36 control subjects were included in the analyses. Relevant subject demographic and injury-related variables are summarized in Table 1. Study groups did not differ in age at MRI, sex, in the number of subjects scanned before vs. after the hardware upgrade, stage of pubertal development, or psychosocial adversity (all P > .05). There were no differences between EI and TBI groups regarding mechanism of injury, time since injury, or ISS. There was a significant difference between EI and TBI groups in ISS excluding head injury (P < .001). As expected, excluding minor injuries, the number of patients with AIS scores ≥2 was higher in the TBI group for head and neck injuries (P < .001) and higher in the EI group for abdominal (P = .009) and extremity (P < .001) injuries.

Global differences

Three-group analyses

Using the MTPC method, there were no significant group X age, group X sex, group X scanner, or group X puberty effects for any graph metrics tested, so interaction terms were removed from the models. Several graph-level metrics were found to differ between groups. Weighted global efficiency (Amtpc = 37,338.1; Acrit = 15,304.3; P < .05) and strength (Amtpc = 72,521.0; Acrit = 14,883.3; P < .05) were significantly lower, and modularity (Amtpc = 37,637.7; Acrit = 2605.6; P < .05) was significantly higher in the TBI group compared to the control group. Strength was also lower in the EI group compared to controls (Amtpc = 31,683.7; Acrit = 12,422.8; P < .05). Modularity (Amtpc = 5606.3; Acrit = 1181.8; P < .05) was significantly higher in the EI group compared to controls, and also higher in the TBI group compared to the EI group (Amtpc = 23,239.7; Acrit = 0; P < .05).

Vertex-level differences

Fig. 1 shows the vertices with significantly higher weighted nodal efficiency for pairwise comparisons across the three groups. Table 2 provides the results of these group comparisons for each region by hemisphere. Although the “control > TBI” comparisons were significant across most regions bilaterally, significant differences were also found for “control > EI” and “EI > TBI” in several (predominantly cortical) regions. Table 3 and Supplementary Fig. 2 (S2) show results of group comparisons for vertex strength. A similar pattern of group differences was noted, with significantly lower vertex strength in TBI than controls in numerous regions within bilateral cortical, cingulate, and subcortical gray matter regions. The EI group showed significantly lower strength relative to controls in bilateral regions, particularly in frontal and parietal lobes. Table 4 and Supplementary Fig. 3 (S3) show regions with significantly lower average path length. Widespread increased path length was found between the TBI group and controls, and was clustered in parietal and temporal lobes in the EI group compared to controls.

Fig. 1

Group differences in vertex weighted nodal efficiency. The vertices shown are those in which weighted nodal efficiency was significantly different between groups. Significance was determined by the MTPC procedure. Vertex color corresponds to lobe membership (red: frontal; green: parietal; blue: temporal; magenta: occipital; yellow: insula; orange: cingulate; light green: subcortical). The left hemisphere is displayed on the left side of each figure.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 2

Multi-threshold permutation correction results for weighted nodal efficiency.

Rows represent brain regions, and columns represent group comparisons. A significant result for the group comparison and region is represented by a”+”, while lack of significant is represented by a”–”.

Region	Control > EI		Control > TBI		EI > TBI
	L	R	L	R	L	R
Frontal
Frontal pole	+	−	+	+	−	+
Lat. orbitofrontal	+	+	+	+	+	−
Med. orbitofrontal	−	−	−	+	−	+
Paracentral	−	−	+	−	−	−
Pars opercularis	−	−	+	−	−	−
Pars orbitalis	−	−	+	+	−	−
Pars triangularis	−	−	−	+	−	−
Rostral middle frontal gyrus	−	−	+	+	−	+
Sup. frontal gyrus	−	−	+	+	−	−
Parietal
Postcentral	−	−	+	−	−	−
Precuneus	−	−	+	+	−	−
Sup. parietal lobule	+	+	+	+	−	−
Temporal
Bank of the sup. Temp. sulc.	−	−	−	+	−	+
Entorhinal	−	−	−	−	+	−
Fusiform	−	+	+	+	−	−
Middle temp. Gyrus	−	−	+	+	−	−
Parahippocampal	−	+	−	+	−	−
Sup. temp. Gyrus	−	−	+	−	−	−
Temp. pole	−	−	−	−	+	−
Occipital
Cuneus	−	+	+	+	−	−
Lat. occipital gyrus	+	+	+	+	−	−
Lingual	−	−	+	+	−	−
Pericalcarine	+	−	+	+	+	+
Insula
Insula	−	−	−	−	+	−
Cingulate
Caudal ant. Cingulate	−	−	+	−	+	−
Isthmus cingulate	−	+	+	+	−	−
Post. cingulate	−	−	+	+	−	−
Rostral ant. Cingulate	−	−	+	+	+	−
SCGM
Amygdala	−	−	−	−	+	−
Caudate	+	+	+	−	−	−
Hippocampus	−	−	+	+	−	−
Nucleus accumbens	−	−	+	−	+	−
Pallidum	−	−	−	−	+	−
Putamen	−	−	+	−	+	−
Thalamus	−	−	+	−	−	−

Table 3

Multi-threshold permutation correction results for vertex strength.

Rows represent brain regions, and columns represent group comparisons. A significant result for the group comparison and region is represented by a”+”, while lack of significant is represented by a”–”.

Region	Control > EI		Control > TBI		EI > TBI
	L	R	L	R	L	R
Frontal
Caudal middle frontal gyrus	−	+	−	+	−	−
Frontal pole	+	−	+	+	+	+
Lat. orbitofrontal	+	+	+	+	+	−
Med. orbitofrontal	+	+	+	+	−	−
Paracentral	+	+	+	−	−	−
Pars opercularis	+	+	+	+	−	−
Pars orbitalis	+	+	+	+	+	−
Pars triangularis	−	+	+	+	−	−
Precentral	+	+	−	+	−	−
Rostral middle frontal gyrus	+	−	+	+	+	+
Sup. frontal gyrus	+	+	+	+	−	−
Parietal
Inf. parietal lobule	+	−	+	−	−	−
Postcentral	+	+	+	+	−	−
Precuneus	−	+	+	+	+	+
Sup. parietal lobule	+	+	+	+	−	−
Temporal
Bank of the sup. temp. sulc.	−	−	−	−	+	+
Entorhinal	+	+	−	−	−	−
Fusiform	−	−	−	+	−	−
Middle temp. gyrus	−	+	+	+	+	+
Sup. temp. gyrus	−	−	+	+	+	−
Transverse temp.	+	+	−	−	−	−
Occipital
Cuneus	−	+	+	+	+	+
Lat. occipital gyrus	+	+	+	+	+	−
Pericalcarine	+	+	+	+	+	+
Insula
Insula	−	+	−	−	+	−
Cingulate
Caudal ant. cingulate	−	−	+	+	+	+
Isthmus cingulate	+	+	+	+	−	−
Post. cingulate	+	−	+	+	+	−
Rostral ant. cingulate	−	+	+	+	+	−
SCGM
Amygdala	+	−	+	−	+	+
Caudate	+	+	+	+	−	−
Hippocampus	−	+	−	+	−	−
Nucleus accumbens	−	−	+	−	+	−
Pallidum	−	+	−	−	+	+
Putamen	−	−	−	−	+	−
Thalamus	+	+	+	+	+	+

Table 4

Multi-threshold permutation correction results for average shortest path length.

Rows represent brain regions, and columns represent group comparisons. A significant result for the group comparison and region is represented by a “+”, while lack of significant is represented by a “–”.

Region	Control > EI		Control > TBI		EI > TBI
	L	R	L	R	L	R
Frontal
Caudal middle frontal gyrus	−	−	−	+	−	−
Frontal pole	−	−	+	+	−	+
Lat. orbitofrontal	+	+	+	+	+	−
Med. orbitofrontal	−	−	+	+	+	+
Paracentral	−	−	+	−	−	−
Pars opercularis	−	−	+	+	−	−
Pars orbitalis	−	−	+	+	+	−
Pars triangularis	−	−	−	+	−	−
Precentral	−	−	−	+	−	−
Rostral middle frontal gyrus	−	−	−	+	+	+
Sup. frontal gyrus	−	−	+	+	−	+
Parietal
Inf. parietal lobule	−	−	−	+	−	−
Postcentral	+	−	+	+	−	−
Precuneus	−	−	+	+	−	−
Sup. parietal lobule	+	+	−	+	−	−
Supramarginal gyrus	+	−	−	−	−	−
Temporal
Bank of the sup. Temp. sulc.	−	−	−	+	−	+
Entorhinal	+	−	+	−	+	−
Fusiform	+	+	+	+	−	−
Middle temp. Gyrus	−	+	+	+	−	−
Parahippocampal	−	+	−	+	−	−
Sup. temp. Gyrus	−	+	−	−	−	−
Temp. pole	+	−	+	−	−	−
Transverse temp.	+	−	−	−	−	−
Occipital
Cuneus	−	+	+	+	−	−
Lat. occipital gyrus	+	+	+	+	−	−
Lingual	−	−	+	+	−	−
Pericalcarine	+	−	+	+	+	+
Insula
Insula	−	−	+	+	−	−
Cingulate
Caudal ant. Cingulate	−	−	+	−	+	+
Isthmus cingulate	−	+	+	+	−	−
Post. cingulate	−	−	+	+	−	−
Rostral ant. Cingulate	−	−	+	+	+	−
SCGM
Amygdala	+	−	+	−	+	−
Caudate	+	−	+	−	−	−
Hippocampus	−	−	+	+	−	−
Nucleus accumbens	−	−	+	−	+	+
Pallidum	−	−	−	−	+	−
Putamen	−	−	+	+	+	−
Thalamus	+	−	+	−	−	−

Network-based statistic

Calculation of the network-based statistic revealed 1 connected component with 22 vertices in which connectivity strength was significantly higher in the control group compared to the TBI group (P = .011). Axial and sagittal views of the vertices and edges in this component are shown in Fig. 2. This component spanned both hemispheres although the majority are left-hemispheric. Furthermore, many of the regions involved were in the temporal and frontal lobes along with subcortical gray matter. Table 5 provides a comprehensive list of the regions in the component with reduced connectivity strength in the TBI group compared to the control group. There were no significant differences between the control and EI group, or between the EI and TBI groups.

Fig. 2

Network-based statistic results testing for group differences.

There was a single connected component with significantly lower connectivity strength in the TBI group compared to the control group. There were no significant differences between the control and EI or the EI and TBI groups.

Table 5

NBS: significant regions for Control > TBI.

Symbols in the table are the same as those in Table 2.

Region	Control > TBI
	L	R
Frontal
Frontal pole	−	+
Lat. Orbitofrontal	+	+
Med. Orbitofrontal	−	+
Pars orbitalis	+	−
Precentral	+	−
Rostral middle frontal gyrus	+	−
Temporal
Fusiform	+	−
Inf. temp. Gyrus	+	−
Middle temp. Gyrus	+	−
Insula
Insula	+	−
Cingulate
Isthmus cingulate	−	+
SCGM
Amygdala	+	−
Caudate	−	+
Hippocampus	−	+
Nucleus accumbens	+	+
Pallidum	+	+
Putamen	+	+
Thalamus	+	−

Associations of graph metrics with age and injury severity

Several graph metrics were associated with age at MRI. Global efficiency (Amtpc = 27,334.0; Acrit = 15,983.0; P < .05) and strength (Amtpc = 26,969.2; Acrit = 17,133.1; P < .05) were positively associated, and characteristic path length (Amtpc = 30,160.4; Acrit = 14,519.7; P < .05) was negatively associated, with age at MRI. Associations with pubertal status (excluding age at MRI from the model) were the same as the results for age at MRI (data not shown).

At the vertex-level, two network metrics were significantly associated with age at MRI in multiple regions. Supplementary Fig. 4 (S4) shows the vertices for which weighted nodal efficiency had a significant positive association with age at MRI (Table 6). Supplementary Fig. 5 (S5) shows the vertices for which average shortest path length had a significant negative association with age at MRI (Table 7).

Table 6

Multi-threshold permutation correction results testing for associations between age at MRI and weighted nodal efficiency.

Symbols in the table are the same as those in Table 2.

Lobe	Region	Pos. effect of age
		L	R
Parietal	Precuneus	+	–
	Sup. parietal lobule	+	+
Temporal	Fusiform	+	–
	Inf. temp. gyrus	+	–
	Middle temp. gyrus	+	+
	Parahippocampal	+	+
Occipital	Cuneus	+	+
	Lat. occipital gyrus	+	–
SCGM	Pallidum	+	+

Table 7

Multi-threshold permutation correction results testing for associations between age at MRI and average shortest path length.

Symbols in the table are the same as those in Table 2.

Lobe	Region	Pos. effect of age
		L	R
Frontal	Pars orbitalis	−	+
Parietal	Inf. parietal lobule	+	+
	Postcentral	−	+
	Precuneus	+	+
	Sup. parietal lobule	+	+
Temporal	Entorhinal	+	−
	Fusiform	+	–
	Inf. temp. Gyrus	+	+
	Middle temp. Gyrus	+	+
	Parahippocampal	+	+
Occipital	Cuneus	+	+
	Lat. occipital gyrus	+	+
Cingulate	Isthmus cingulate	+	−
	Post. cingulate	+	−
SCGM	Pallidum	+	+
	Thalamus	−	+

We also tested for associations between graph-level metrics and injury severity, both the ISS and lowest GCS scores. There were no significant correlations with ISS for any network metrics, nor were there any significant group X ISS interactions. Global efficiency (Amtpc = 19,984.0; Acrit = 16,072.9; P < .05) and strength (Amtpc = 19,728.9; Acrit = 16,908.0; P < .05) were positively associated with lowest GCS.

Discussion

Despite the high incidence of pediatric TI, little is known regarding how different types of injury affect overall brain network organization and its individual components. Here, we used graph theory methods to analyze DTI tractography networks in children approximately 2 months post-TI, and found widespread differences between children with TBI, EI, and healthy controls without injury. The effect of group on all network metrics was consistent across age, pubertal stage, and sex. At the graph (whole brain) level, we found that strength was lower and modularity was higher in both the EI and TBI groups compared to the control group suggesting more segregated networks following TI. Furthermore, global efficiency was lower in the TBI compared to the control group, which is indicative of a less integrated network. At the vertex (regional) level, connectivity strength and nodal efficiency were significantly lower for the TBI group compared to both EI and control groups, and lower in the EI compared to the control group. Lower connectivity strength and nodal efficiency was widespread such that TI effects were evident in brain regions spanning all major lobes of the brain and in both hemispheres. EI-TBI differences were particularly prevalent in prefrontal-limbic circuitry including bilateral frontal, subcortical/striatal, and cingulate regions. Vertex-level average shortest path length similarly showed extensive group differences, but was significantly higher in TBI compared to the EI and control groups. Despite the extensive significant group differences in vertex-level graph metrics, the NBS analysis did not show any differences between the EI group and either the TBI or control group. There was a single connected component with significantly lower connectivity in the TBI group compared to controls, largely consisting of left-lateralized frontal and subcortical regions. These results together show that there are quantifiable differences between non-injured children and children with TI, as well as between children with TBI and EI, which are in agreement with our hypothesized relationships of injury type to network structure. Increasing age was associated with increases in global efficiency and strength and decreases in characteristic path length. The increased segregation associated with TI suggests that injured children may have greater reliance on neural subnetworks and reduced interconnectivity with other regions, which may contribute to the cognitive and psychological health difficulties often experienced after TI in general and TBI in particular.

Graph theory analysis is well-suited to study brain structural changes after pediatric TI and has been applied to studies of TBI but not EI. Consistent with our finding of decreased global efficiency and increased modularity in TBI versus controls, Yuan et al., 2015 found a similar pattern in children with mild TBI (mTBI) within the first week following injury (Yuan et al., 2015). A combination of higher efficiency and lower modularity, as seen in the control group, has been shown to provide an ideal network structure for performing more cognitively-demanding tasks (e.g., working memory) (Kitzbichler et al., 2011). Our findings suggest that this disruption to the brain structural network following TBI persists beyond the acute stage of recovery, and may contribute to the deficits in working memory and complex cognitive processes present following TBI (Babikian et al., 2011). To our knowledge, only one other group has reported regional network metrics across the whole brain, showing lower local efficiency predominantly in frontal and occipital regions, similar to our results comparing nodal efficiency (Caeyenberghs et al., 2012). Yuan and colleagues analyzed only a subset of brain regions but found lower degree (number of connections) in several frontal and occipital regions, which is analogous to the differences we found in regional strength (Yuan et al., 2015).

Effects of trauma on brain connectivity

While other studies of TI have not investigated structural brain network differences between participants with brain and bodily injury in relation to healthy controls, several studies have reported the effect of traumatic exposures (related to bodily threat from natural disasters or from maltreatment) on network metrics. Suo and colleagues examined children exposed to a natural disaster who developed PTSD in relation to trauma-exposed children who did not (Suo et al., 2017). Graph theory analysis of DTI data revealed that the PTSD group had lower global efficiency and increased path length compared to controls, indicating reduced network integration. Moreover, in children with previous exposure to maltreatment, vertex strength was significantly lower throughout the frontal and temporal lobes compared to controls (Puetz et al., 2017). The impact of childhood stress exposure on WM appears to persist into adulthood, as well. In young adults with moderate-to-high prior exposure to maltreatment, strength and global efficiency were found to be lower compared to young adults with no-to-low levels of exposure (Ohashi et al., 2017). Taken together, these studies converge on the identification of concurrent increase in segregation and reduction of integration of brain sub-networks in children exposed to a wide range of traumatic stressors compared to healthy children without trauma or injury. In combination with our findings, these results suggest that traumatic stress exerts a significant impact on the structural characteristics and connectivity of widespread sub-networks that cannot be explained solely by direct trauma to the brain. It is possible that the activation of HPA and/or noradrenergic stress systems identified following TI (Ewing-Cobbs et al., 2017) contributes to the alteration in brain network architecture following TI. When testing for network differences associated with severity within the TBI group only, our findings support the hypothesis that more severe TBI (as measured by the lowest GCS score) leads to larger changes in the whole-brain network.

Developmental impact of TBI

Normal human brain development is marked by changes in gray matter volume and thickness and white matter connectivity which proceeds along a posterior-to-anterior axis, such that frontal regions and WM tracts continue to undergo changes through adolescence (Lebel and Beaulieu, 2011; Sowell et al., 2004). At the global network level, efficiency increases and modularity and characteristic path length decrease over time (Cao et al., 2014; Dennis et al., 2013; Hagmann et al., 2010; Wierenga et al., 2016). At the regional level, connectivity strength increases in frontal regions into late adolescence (S. T. Baker et al., 2015). The network-based effects identified in the present study tend to be aggregated in later-developing fronto-temporal regions, consistent with research showing increased vulnerability to disruption in brain regions undergoing rapid development (Babikian et al., 2015). These results also resemble those of our prior study, based on a different sample of children and adolescents with TI, which highlighted the vulnerability of association pathways (particularly those in the left hemisphere) that connect prefrontal regions with short- and long-range fibers from other lobes (Ewing-Cobbs et al., 2016). The greater disruption of connectivity in the left than right hemisphere specifically following TI may be related to the protracted development of left hemisphere association pathways during childhood and adolescence (Bonekamp et al., 2007; Simmonds et al., 2014; Wilde et al., 2010; Yu et al., 2014). Additional work is needed to untangle the effects of injury at different stages of chronological and pubertal development on regional and global connectivity metrics.

Limitations

Our study has some limitations. The scanner was upgraded in the middle of the study; we included a covariate in our statistical models to account for potential differences due to this change. Additionally, this is only a cross-sectional analysis and we cannot comment on developmental changes in brain network metrics in children with TI. Our sample was injured in motor vehicular incidents and received treatment or hospitalization. Therefore, our findings may not generalize to children with different types of injury mechanisms or those not seeking treatment. Furthermore, subjects with large lesions were excluded from analyses, which may underestimate the true effects of TBI on brain network connectivity since those with the largest lesions may have greater disruption of connectivity. We were conservative in our quality assessment for inclusion in the study, leading to the exclusion of many subjects due to movement-related artifact. A higher proportion of TBI subjects was excluded for this reason; this reduction in overall sample size reduces power to identify true group differences. Based on the heterogeneity in TI patients, specifically the unequal distribution in injury severity subgroups, a larger sample would provide greater power to identify differences among these groups. It is possible that the findings in some of the TI participants reflect preinjury differences in organization of brain networks rather than post-injury changes.

Conclusions

In summary, we found widespread differences in brain network metrics consistent with lower integration and higher segregation in children with TI compared to a control group. These findings parallel the pattern of deficits in cognitive and behavioral outcomes (e.g., memory deficits, emotional regulation) due to both TBI and non-head injury reported in large, longitudinal studies of children with TI (Babikian et al., 2011; Fay et al., 2009). The consistent deficits in structural connectivity among and between fronto-limbic structures may contribute to the cognitive and behavioral changes found following TI. Taken together, we show that TI results in a global disruption of structural brain networks which is more extensive with more severe injury. Future research in this cohort will investigate functional brain network structure, will assess network changes over the course of 1 year, and will explore associations between changes in brain network metrics and cognitive and behavioral outcomes. Of particular interest will be possible divergence in network organization between children with different types of injury and levels of injury severity, and whether such changes can predict long-term neurobehavioral outcomes.

Conflicts of interest

The authors report no financial or other conflict of interest.

Funding

This work was funded by the National Institutes of Health grant R01-NS046308. The content is solely the responsibility of the authors and does not necessarily represent the official views of the granting institute. Support was also provided by the Nicole and Evan Katz Pediatric Neurodevelopmental Fund.