Network topology changes in chronic mild traumatic brain injury (mTBI).

BACKGROUND: In mild traumatic brain injury (mTBI), diffuse axonal injury results in disruption of functional networks in the brain and is thought to be a major contributor to cognitive dysfunction even years after trauma.
OBJECTIVE: Few studies have assessed longitudinal changes in network topology in chronic mTBI. We utilized a graph theoretical approach to investigate alterations in global network topology based on resting-state functional connectivity in veterans with chronic mTBI.
METHODS: 50 veterans with chronic mTBI (mean of 20.7 yrs. from trauma) and 40 age-matched controls underwent two functional magnetic resonance imaging scans 18 months apart. Graph theory analysis was used to quantify network topology measures (density, clustering coefficient, global efficiency, and modularity). Hierarchical linear mixed models were used to examine longitudinal change in network topology.
RESULTS: With all network measures, we found a significant group × time interaction. At baseline, brain networks of individuals with mTBI were less clustered (p = 0.03) and more modular (p = 0.02) than those of HC. Over time, the mTBI networks became more densely connected (p = 0.002), with increased clustering (p = 0.001) and reduced modularity (p < 0.001). Network topology did not change across time in HC.
CONCLUSION: These findings demonstrate that brain networks of individuals with mTBI remain plastic decades after injury and undergo significant changes in network topology even at the later phase of the disease.

Introduction

Traumatic Brain Injury (TBI) is a debilitating neurological disorder that results from a traumatic external impact to the head or body. TBI is a significant public health burden, affecting 10 million people worldwide every year (Humphreys et al., 2013, Hyder et al., 2007). Veterans and those who serve in the military are at an especially high risk. Data from the Department of Defense revealed that from 2000 to 2011, approximately 4.2% of those who served in the armed forces were diagnosed with some form of TBI (Centers for Disease Control and Prevention, 2015). TBI severity is typically defined as mild, moderate, or severe, based on time spent unconscious and/or coma rating score, the duration of post-traumatic amnesia, and neuroimaging results. Of these categorizations, by far the most common is mild traumatic brain injury (mTBI), which accounts for 70–90% of all TBI cases (Leo et al., 2015).

Though on the lower end of the severity spectrum, mTBI often produces a wide range of neurological impairments that can have long-lasting consequences on brain health (Dean and Sterr, 2013, Quinn et al., 2018). These impairments include cognitive disruptions, motor control deficits, and behavioral problems (Basford, 2003, Rabinowitz and Levin, 2014). Approximately 20–30% of mTBI cases (Shenton, 2012) and up to 65% of moderate/severe cases (Rabinowitz and Levin, 2014, Selassie, 2008) report chronic problems that can be detected years after initial recovery. These persistent deficits cause disability and interfere with a patient’s ability to perform day-to-day tasks (Rabinowitz and Levin, 2014). To date, very few studies have investigated how brain networks adapt to the chronic progression of mTBI (Caeyenberghs et al., 2017), resulting in a limited understanding of how the brain reacts to this type of insult over a long period of time. To better elucidate this question, we conducted a longitudinal study of chronic mTBI to assess changes in brain network topology over an 18-month period.

Functional magnetic resonance imaging (fMRI) allows the examination of neuronal activity by tracking its metabolic and hemodynamics consequences (Fox and Raichle, 2007), and has thus provided critical insights into brain function (Biswal et al., 2010, Nakamura et al., 2009, Raichle, 2011). With respect to mTBI, resting state fMRI (rsfMRI) studies have shown abnormal patterns of activity within the default-mode network (DMN) (Mayer et al., 2015, Zhou et al., 2012), as well as other resting-state networks, including the motor and fronto-parietal networks (Palacios et al., 2017, Stevens, 2012, Zhou, 2014). In recent years, functional connectivity alterations associated with TBI have been examined using the mathematical concepts of graph theory (Caeyenberghs et al., 2017). Graph theory provides a powerful framework for characterizing the topology of complex biological systems (Bassett and Bullmore, 2009, Sporns, 2018, Sporns, 2013). In this perspective, brain networks can be thought of as graphs which are composed of sets of nodes (e.g. distinct brain regions), linked by sets of edges, which can be either structural (white matter fiber tracts) or functional (correlated activity between regions) in nature. Graph analysis thus allows an integrative approach to characterizing different properties of brain networks, including integration, segregation, and node centrality.

Recent studies have shown graph theory to be a promising method to model brain network organization and to quantify pathological deviations across a variety of disorders (Caeyenberghs et al., 2017, Ajilore, 2014, Bassett et al., 2012). In the case of TBI, there have been several imaging studies that have used this approach to characterize alterations in brain network topology in response to trauma (for reviews see (Caeyenberghs et al., 2017, Imms, 2019). These studies indicate that TBI is associated with network hyperconnectivity, characterized by increased density and clustering coefficient, and suboptimal global integration (Caeyenberghs et al., 2017, Aerts et al., 2016). While these are useful insights, it is important to note that there is a high degree of heterogeneity in findings (Caeyenberghs et al., 2017). It is likely that certain aspects which are unique to TBI research, such as difference in severity and time from injury, contribute to the divergent nature of findings across studies. It thus becomes important to distinguish between studies that have looked at mild vs. severe/moderate TBI, as well as those that have studied acute vs. chronic progression of the illness.

One specific subgroup of TBI that has received relatively limited attention is chronic mTBI. While TBI is a static event, its etiology is better thought of as a chronic disease process (Masel and DeWitt, 2010). There is evidence that trauma triggers persistent neuroinflammation and long-term atrophy of both grey and white matter, leading to a progressive degeneration of the neural substrate (Bigler, 2013, Farbota, 2012, Gilmore et al., 2020). For example, in a longitudinal study, Farbota (2012) followed individuals for up to 4 years post injury and found that, rather than showing just an acute period of brain degeneration, TBI resulted in a protracted period of brain change that lasted for years following trauma. These findings raise important questions about the chronic consequences of mTBI and indicate that neural degeneration is a long-term process with a time course that has not yet been clearly elucidated.

To our knowledge, only two longitudinal studies using graph analysis of rsfMRI have been conducted in mTBI. Both studies used rsfMRI to assess brain networks immediately in the acute phase of injury and again at 6 months and 1 year respectively (Dall’Acqua et al., 2017, Messé et al., 2013). Messé et al. (2013) found significantly lower network modularity in patients who suffered from post-concussive symptoms compared to those that did not. No other differences in network properties were identified. Dall’Acqua et al. (2017) found that the acute phase of recovery was characterized by functional hypoconnectivity in the DMN (reduced network strength between controls and mTBI patients), and that this hypoconnectivity normalized over the course of a year.

These studies offer crucial, if limited, insights into the restructuring of brain networks after mTBI. The present study expands on these findings by investigating changes in network topology in chronic mTBI (~20 years post trauma on average). This is the first longitudinal study using graph theory-based analysis to investigate the progression of mTBI at such a late phase of injury. The results provide an insight into how brain network topology continues to change even decades after initial recovery.

Methods

Participant demographic and clinical information

Participants were veterans receiving services in the Minneapolis VA Health Care System (MVAHCS). Participants in the TBI group were diagnosed with mild or symptomatic TBI based on the Mayo TBI Severity Classification System (Malec et al., 2007) following information self-reported during an interview using the Minnesota Blast Exposure Screening Tool (MN-BEST) (Nelson et al., 2015). Participants in the healthy control group (HC) were veterans with no history of TBI or head trauma. Informed consent was obtained after the study was approved by the MVAHCS Institutional Review Board. A total of 139 veterans were enrolled in the larger longitudinal study (see (Gilmore et al., 2020) for overall study and sample details). Of these, 50 in the mTBI group and 40 in the HC group were able to undergo functional MRI testing on two occasions separated by 18 months. Participants completed a clinical interview at their baseline visit that included the Mini-International Neuropsychiatric Interview (M.I.N.I.)(Amorim, 2000), Alcohol Use Disorders Identification Test (AUDIT-C)(Reinert and Allen, 2002) and the MN-BEST.

Mild TBI severity scoring

Ratings of mTBI likelihood and severity were assigned by consensus of doctoral level clinical neuropsychologists based on information secured by trained study interviewers using the semi-structured MN-BEST. Symptoms of mTBI assessed included altered consciousness (e.g., confusion and disorientation), loss of consciousness (LOC) <30 min, post-traumatic amnesia (PTA) up to 24 h, and neurological symptoms (e.g., headache, tinnitus, nausea, sensitivity to light or noise) immediately after the event. The three most significant potential blast-related and impact head injury events were considered, each of which received a severity score ranging from 0 (no concussion) to a potential maximum of 30 (severe TBI). No score was higher than 4 (the maximum within the mTBI range) for a single event in the current sample. The severity scores for each of the reported head injury events were summed to get an overall mTBI severity score for each participant.

Neuropsychological assessment

A battery of standard tests were used to examine neuropsychological symptoms known to be impacted by TBI. The Depression Anxiety and Stress Scale (DASS-21) (Ownsworth et al., 2008) was used to assess affect. The Neurobehavioral Symptom Inventory (NSI) examined post-concussion symptoms(Belanger, 2017), and the Post-Traumatic Stress Disorder checklist (PCL) (Blevins et al., 2015) was used to examine PTSD symptoms. All scales were completed at baseline and at the follow up visit prior to the second scan.

Image acquisition and preprocessing

MRI data were acquired at the University of Minnesota’s Center for Magnetic Resonance Research on a single 3 T Siemens Prisma fit scanner (Siemens, Erlangen, Germany) equipped with 32-channel head coils. For each participant, T1-weighted (MP-RAGE, TR = 2400 ms, TE = 2.24 ms, TI = 1060 ms, voxel size = 0.8 mm isotropic, flip angle = 8 degrees) and T2-weighted (SPACE sequence, TR = 3200 ms, TE = 564 ms, voxel size = 0.8 mm isotropic, variable flip angle) scans were acquired. A resting state fMRI (rsfMRI) scan was acquired after the participants were instructed to keep their eyes open, look at a fixation cross, and not think about anything in particular. The Human Connectome Project (HCP) aging and development protocol (HCP-A/D, Harms et al., 2018) (TR 800 ms, TE = 37 ms, voxel size 2 mm isotropic, multiband acceleration factor = 8, 600 volumes, 8 min) was used for acquisition.

Image quality for the rsfMRI data was using tools from the FSL toolbox (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/; v6.0.1) and custom MATLAB (R2017b) scripts. Preprocessing for the image quality assessment consisted of motion correction and spatial smoothing (6 mm full width at half maximum). Framewise displacement (FD) and root mean variance of temporal derivative over voxels (DVARS) was calculated using the methods outlined in Power et al. (Power et al., 2012). Volumes with FD > 0.5 mm and/or DVARS > 8 were considered to have suspect quality, along with 1 vol before and 2 volumes which followed. The total number of flagged volumes within a rsfMRI scan was then computed; runs with>25% flagged volumes were rejected from further analysis.

Following acquisition, the scans were pre-processed using a container of version 4.0.1 of the HCP minimal preprocessing pipeline (Glasser et al., 2013). The container can be made available upon request. T1-weighted and T2-weighted anatomical scans were aligned, and gradient distortion correction was applied. Segmentation and cortical reconstruction of the structural images to the cortical sheet was completed with FreeSurfer (v6.0.0) (Dale et al., 1999, Fischl et al., 2001, Fischl et al., 1999). Structural images were then output to CIFTI space (Glasser et al., 2013). Preprocessing of the rsfMRI scans involved correcting for nonlinear gradient inhomogeneities and realignment of EPI volumes to correct for participant head motion by registration to a single band reference image (using the FSL tool “FLIRT” Jenkinson et al., 2002, Jenkinson and Smith, 2001). To correct for distortion in the phase encoding direction, a pair of spin echo images with AP and PA phase encoding polarity (with matching geometry and echo-spacing to the rest scan) was also acquired. These spin echo images were used to estimate a distortion field with the FSL tool “topup” (Andersson et al., 2003). This estimated distortion field was then applied to the rsfMRI volumes with “FLIRT”. The rsfMRI images were then registered to the T1-weighted and T2-weighted images, and finally to MNI space. To complete the minimal pre-process of the rsfMRI scan, its MNI registered volumes were resampled onto cortical surfaces and extracted to CIFTI space.

The final major step, upon completion of both the structural pipeline and minimal functional pipeline, was to correct temporal artefacts within the rsfMRI time series. First, the rsfMRI time series was linearly detrended with a high band pass filter cutoff of 2000s with a slow roll off (Smith et al., 2013). The fMRI data, including the rs-fMRI, was then concatenated across scans within a given scanning session so as to provide more data as input to perform spatial independent component analysis with the FSL tool “MELODIC” and “FIX” (Griffanti et al., 2014, Salimi-Khorshidi et al., 2014) and to yield better separation of “signal” and “noise” components. Once “FIX” had classified the spatial components as “signal” or “noise”, both component types were regressed into the data and the resulting “noise” spatial maps were multiplied by the associated time series and subtracted from the original dataset. Finally, following the ICA + FIX process, surface-based functional alignment was run on the ICA + FIX cleaned data with the tool MSMAll (Robinson et al., 2014). MSMAll employed myelin maps, rsfMRI network maps, and rsfMRI visuotopic maps to align a participant’s cortical data to a group template. Once complete the MSMAll and ICA + FIX fMRI concatenated data was dissociated into each respective fMRI scan and this rsfMRI data was used for subsequent analyses. Prior to the final analysis, individuals who exhibited > 25% flagged volumes (image volumes deemed to have an excess amount of noise artifacts) were dropped from the sample. The final sample size was 42 mTBI and 33 HC. We found no significant group differences in the number of subjects dropped or percent noisy volumes.

Brain network construction

Graph theory analysis was applied to resting state fMRI data using a multivariate approach (Fig. 1). Functional connectivity was computed by taking the absolute value of the Pearson's correlation between all possible pairs of time series, creating a 419 × 419 (N × N) connectivity matrix based on a combination of the 400 ROI Schaefer cortical (Schaefer et al., 2018) and 19 ROI Harvard-Oxford subcortical (Desikan et al., 2006, Frazier, 2005, Makris, 2006) atlases using the tool MACCHIATO (https://github.com/tjhendrickson/MACCHIATO). The raw functional connectivity matrices for each subject were then thresholded and binarized to facilitate network topology analysis. We applied a thresholding method for each graph which kept only statistically significant connections between nodes. This was achieved by testing the null hypothesis that the correlation coefficient at each link was significantly different from zero between any node pairs. Two nodes were considered functionally connected if their correlation coefficient was significant at p < 0.05 level, resulting in an adjacency matrix G containing information about the connectivity structure of the network. Correction for multiple comparisons over all entries in the matrices for each session was performed using the linear step-up procedure for false discovery rate (FDR). After thresholding, networks were then binarized by assigning a 1 to significant links (edges) and a 0 to non-significant links.

Fig. 1

Schematic representation of graph theoretical analysis using resting-state functional MRI data. After data acquisition (A) and preprocessing (B), the fMRI data was parcellated (C) and time course data was extracted from each region of interest (D). A 419 × 419 matrix containing functional connectivity data was computed by taking the absolute value of the Pearson’s correlation between all possible pairs of time series (E). To reduce complexity and generate sparse graphs for network analysis, a threshold was applied to the functional connectivity matrix which kept only significant links (connections) between nodes. Links which survived thresholding were assigned a value of 1, and those that did not were assigned 0′s to generate unweighted binarized networks (F). Binarized networks were used to compute network topology metrics (H).

Graph theory analysis

The Brain Connectivity Toolbox (Rubinov and Sporns, 2010) and our own custom scripts were used to compute graph metrics on binarized networks. Networks were analyzed using multivariate graph metrics which characterize global network topological features such as global integration and functional segregation. Network density and global efficiency were used to assess global integration, while clustering coefficient and modularity were used to assess functional segregation. The graph theory measures investigated in the study considered a network , composed of a set of nodes and a set of edges (significant connections between node pairs).

Network density refers to the fraction of total connections in a network compared to the total number of possible unique connections(Rubinov and Sporns, 2010):

Global network efficiency is based on characteristic path length (defined as number of steps, or links, between network nodes) and measures how efficiently a network can exchange information (Latora and Marchiori, 2001). The global efficiency, , is defined as the average nodal efficiency, where corresponds to the path length between nodes and :

Clustering Coefficient (C represents the degree of local connectivity of a given node with its nearest neighbors ():

Here, represents the degree of node and represents the number of edges between the neighbors of node The average clustering coefficient for the entire network is the average clustering coefficient at each node (Watts and Strogatz, 1998):

Modularity characterizes the degree to which nodes are partitioned into a set of modules or communities having dense interconnected nodes and less communication with nodes of the other modules (Newman, 2006). Here, we use the Louvain method for community detection (Blondel et al., 2008) to find the optimal partition that maximizes the modularity value. The measure of modularity is defined as:where is the number of connections attached to the node (or node degree), is the partition the node belongs to and and otherwise.

Statistical analysis

Statistical analysis was performed in R (v3.6.2). Neurobehavioral and clinical measures (NSI, PCL, DASS-21) were assessed using general linear models, testing for the main effects of group, time and a time × group interaction. Graph metrics were analyzed using a hierarchical linear modeling approach using the nlme package (Pinheiro et al., 2020). In the final model, we included the main effects of group (TBI vs. HC), time (baseline vs. follow-up), and the interaction between group and time. We also included age (age at baseline scan) and the number of flagged volumes (for image quality) in the final model. Subject-specific random intercepts were included to account for repeated measures across time for each subject within each group. The model selection process involved the inclusion of additional data features in the initial model, including, gender, time since injury, TBI severity, PTSD symptoms, and level of education. However since these factors were not significant and did not improve model fit (based on Akaike information criterion Akaike, 1974), they were excluded from the final model. Normality of model residuals was confirmed by visual inspection of residual distribution (Casson and Farmer, 2014) and the shapiro-wilks test for normality (Shapiro et al., 1965) (p = 0.13).

Planned contrasts were used to further explore significant main effects and interactions. Contrasts were computed using the least-square mean approach for each level of group and time, resulting in 6 orthogonal post-hoc tests; FDR correction was used to adjust for multiple comparisons in all tests.

Results

Clinical demographics & measures

A total of 42 veterans in the mTBI group and 33 in the HC group were included in the final analysis. The mTBI and HC groups were similar with regard to age, gender, education level, and prevalence of hypertension, diabetes, hazardous drinking/alcohol use disorder, and PTSD (Table 1). The average time from injury in the mTBI group was 20.7 years, with a large standard deviation of 15.6 years. As expected, we observed a significant difference across groups for neurobehavioral symptomatology (χ(Hyder et al., 2007) = 5.94, p = 0.011), with the mTBI group showing worse symptoms on the NSI total score. There was no change across time in either group. There was a trend level group effect on the PCL total score (χ2 = 3.41, p = 0.063), with increased scores in the mTBI group compared to controls. PCL total score did not change as a function of time (χ2 = 0.083, p = 0.773). Similarly, we found no significant differences between groups on the DASS-21 subscales measuring stress, anxiety, and depression (p > 0.05 for all comparisons).

Table 1

Participant characteristics.

	Healthy controls	mTBI	Statistical test
N	33	42	–
Age: mean (SD)	47.0 (11.9)	50.9 (12.1)	t = −1.59, p = 0.11†
Women: N (%)	9 (27%)	5 (12%)	χ2 = 2.87, p = 0.09
Mean Education Levela: mean (SD)	6.5 (1.5)	6.7 (1.8)	χ2 = 0.22, p = 0.78
Diabetes: N (%)	7 (21%)	8 (19%)	χ2 = 0.17, p = 0.83
Hazardous Drinking/AUDb: N (%)	13 (39%)	15 (36%)	χ2 = 0.11, p = 0.74
PTSDc: N (%)	1 (3%)	3 (9%)	χ2 = 0.61, p = 0.43
Years from TBI: mean (SD)	–	20.7 (15.5)	–
MN-BEST Severity Total Score: mean (SD)	–	1.88 (2.52)	–

TBI: Traumatic Brain Injury; AUD: Alcohol Use Disorder; PTSD: Post-Traumatic Stress Disorder

a. Education level was coded as follows: 1 = <8 years, 2 = Some High School, 3 = Graduated High School, 4 = GED, 5 = Work towards Associate's Degree, 6 = Completed Associate's Degree, 7 = Work towards Bachelor's Degree, 8 = Completed Bachelor's Degree or higher, 9 = Other.

b. Rates of possible hazardous drinking/AUD were determined by the AUDIT-C (a score of ≥ 3 for women and ≥ 4 for men).

c. Rates of PTSD were determined by meeting criteria for current (within the past month) PTSD on the M.I.N.I.

Statistical test used were either students t-test or a chi-square test.

Global integration

With respect to network density (the fraction of present connections to all possible connections), we observed a significant effect of group (χ2 = 5.16, p = 0.022) and a significant group × time interaction (χ2 = 7.15, p = 0.007; Fig. 2; Supplementary Table 2). At the baseline scan, there was a trend towards a group difference in density (t = 2.31, p = 0.064), with a sparser network topology in the mTBI group. Within the mTBI group, we found a significant difference across time (t = −3.82, p = 0.002; d = 0.630), with density increasing over time. There was no significant change in density in HC (t = 0.152, p = 0.879). We also found a significant main effect of age on network density (χ2 = 13.2, p > 0.001).

Fig. 2

Network density represents the fraction of total connections in a network compared to the total number of possible unique connections. Average network density across the two timepoints is plotted for healthy controls (left panel) and for individuals with mTBI (right panel). Asterisks represent a significant time effect within groups.

With global efficiency we found a significant interaction between group and time (χ2 = 6.95, p = 0.008), with global efficiency increasing over time in the mTBI group (t = −3.75, p = 0.001; d = 0.590; Fig. 3, Supplementary Table 3). No change in global efficiency was seen in HC. A trend level difference in global efficiency between groups was observed at baseline scan (t = 2.51, p = 0.061). Interestingly, we also found a significant effect of age (χ2 = 13.7, p < 0.001), indicating that higher age was associated with reduced global efficiency.

Fig. 3

Global network efficiency is based on characteristic path length (defined as number of steps, or links, between network nodes) and measures how efficiently a network can exchange information. Global efficiency across the two timepoints is plotted for healthy controls (left panel) and for individuals with mTBI (right panel). Asterisks represent a significant time effect within groups.

Functional segregation

For clustering coefficient, we found a significant group effect (χ2 = 6.18, p = 0.012) and a significant group × time interaction (χ2 = 6.77, p = 0.009; Fig. 4). At baseline, the mTBI group showed reduced clustering compared to HC (t = 2.48, p = 0.030; d = 0.674). Over time however, clustering significantly increased in the mTBI group (t = −3.90, p = 0.001; d = 0.665). There was no change in clustering in the brain networks of HC (t = −0.015, p = 0.988).

Fig. 4

The average clustering coefficient quantifies the abundance of connected triangles in a network and characterized how densely connected nodes are to their immediate neighbors. Average clustering coefficient across the two timepoints is plotted for healthy controls (left panel) and for individuals with mTBI (right panel). Asterisks represent a significant time effect within groups.

With modularity, we found a significant group effect (χ2 = 6.44, p = 0.011), and a significant group × time interaction (χ2 = 8.26, p = 0.004; Fig. 5). Contrasts revealed a significant between group difference in modularity at baseline (t = −2.53, p = 0.027; d = 0.650), with the mTBI group showing elevated modularity compared to HC. In the mTBI group modularity decreased over time (t = 4.28, p < 0.001; d = 0.704), with no change in HC (t = −0.015, p = 0.989).

Fig. 5

Modularity characterizes the degree to which nodes are partitioned into a set of modules or communities which have dense interconnected nodes and less communication with nodes of the other modules. Modularity across the two timepoints is plotted for healthy controls (left panel) and for individuals with mTBI (right panel). Asterisks represent a significant time effect within groups.

For both clustering coefficient and modularity, we found a significant effect of age. For clustering coefficient (χ2 = 7.51, p = 0.006), increased age was associated with decreased clustering, whereas with modularity (χ2 = 11.7, p < 0.001), increased age correlated with increased modular network structure.

For all graph metrics, we found a significant main effect of the number of flagged volumes (see Supplementary tables). There was no difference in the number of flagged volumes across groups (χ2 = 0.176, p = 0.674), or across time (χ2 = 0.919, p = 0.337).

Discussion

We utilized graph theory analysis to examine how resting state functional connectivity would change over time in a large sample of 42 veterans with mTBI and 33 controls. To the best of our knowledge, this is the first study to longitudinally track changes in brain network topology in chronic mTBI. The unique aspect of our sample was that individuals with mTBI were, on average, 20.7 years removed from the time of injury. No studies have examined brain network changes in mTBI at such chronic timescales. Our results demonstrate three main findings: i) brain networks of individuals with mTBI remain plastic decades after injury and undergo significant changes in network topology, ii) global integration increased over time in mTBI group to the level of HC, and iii) brain networks of individuals with mTBI became more clustered and less segregated into modules over time.

Global integration metrics, such as global efficiency and network density, characterize the brain’s ability to rapidly combine specialized information from distributed brain regions, and disruptions in these metrics have been observed in both functional and structural networks in TBI (Caeyenberghs et al., 2017, Caeyenberghs, et al., 2012, Han et al., 2016). Alterations in global integration have been tied to disruptions in long range fibers that can result from traumatic axonal injury as demonstrated in diffusion studies (Kuceyeski et al., 2016, Pandit et al., 2013). For example, Pandit et al. (2013) found that individuals with TBI who showed greater levels of white matter damage (specifically to long-range connections), exhibited longer characteristic path lengths as well as reductions in network density and global efficiency. Previous longitudinal studies examining global integration in TBI have reported diverging results. In a small study of severe TBI (n = 6), global integration metrics (global efficiency and path length) were higher in TBI compared to controls at 3 months post injury but decreased significantly over the course of 6 months (Nakamura et al., 2009). In a larger study of mTBI (n = 17), no differences between healthy controls and mTBI were found in respect to global integration (Messé et al., 2013).

Given the differences in population and late phase of injury, it becomes difficult to compare these findings to our own results. In our sample, when assessing global efficiency, we found a trend level difference at baseline (p = 0.06), with the mTBI group showing reduced efficiency compared to HC. Over time, we saw a significant positive shift in global efficiency in the mTBI group, with no change in HC. A similar pattern emerged with network density, where individuals with mTBI showed lower network density at baseline but a significant increase over time. The increase in overall connectivity could help explain why global efficiency increased in the mTBI group, since more connections (i.e., higher density) might result in more routing paths and possibly more efficient information propagation. Because global efficiency mainly reflects integrative information processing across distal brain regions (Sporns, 2018, Ajilore, 2014, Desikan et al., 2006), the observed increase in efficiency may imply an adaptation that results in a higher degree of long range connections in mTBI (Bassett and Bullmore, 2009). Though this pattern of change serves to “normalize” levels of global integration to the level of controls, increased integration should not automatically be seen as supporting positive compensation in brain injury. Rather, it can be hypothesized that increases in long range connections are likely metabolically costly and thus potentially maladaptive (Casson and Farmer, 2014, Shapiro et al., 1965, Caeyenberghs, et al., 2012). Indeed, this line of thinking is supported by research that associates functional hyperconnectivity with reduced cognitive performance in TBI (Bassett and Bullmore, 2009, Newman, 2006).

The increased global efficiency of functional brain networks in mTBI may not be entirely due to more connections, but may also reflect an intrinsic alteration in brain wiring patterns. Examining changes in functional segregation may provide insights to how these alterations emerged. Functional segregation refers to the ability of networks to segregate neuronal processing amongst functionally related regions that are arranged in modules (Sporns, 2018, Sporns, 2013). The clustering coefficient is a common measure of functional segregation and serves to characterize the density of a node's connections to its neighbors. At the baseline scan, we found a significant group difference in clustering, with the mTBI group showing reduced clustering compared to HC. Over time however, clustering increased in the mTBI group to the level of controls (no significant time effect in HC). Several studies have shown increased clustering in TBI compared to controls (Bassett and Bullmore, 2009, Newman, 2006, Han et al., 2016), and these findings have been interpreted to support a ‘hyperconnectivity’ hypotheses of recovery from TBI (Caeyenberghs et al., 2017). However, still other studies have failed to find differences in clustering (Nakamura et al., 2009, Pinheiro et al., 2020), challenging the validity of this claim. Our findings do not directly support a hyperconnectivity hypothesis since clustering increased only to the level of controls.

Modularity is another useful tool to characterize network topology. Modularity describes the propensity of a network to subdivide into modules that are defined by a high density of intra vs. inter modular connections. Individuals with mTBI had higher modularity compared to controls at baseline, but underwent a significant negative shift over time. These findings suggest that at baseline, the brain networks in mTBI were less interconnected (lower global integration, lower clustering) and more subdivided into modules compared to HC. Over time, however, there was a shift in topology to a more interconnected and less modular structure. The shift is especially striking given the lack of change in the networks of HC.

The negative shift in modularity is consistent with the positive shift in global efficiency, as less modular networks tend to have more long-range connections that subserve global information transfer (Bullmore and Sporns, 2012, Di et al., 2013). The economic theory of brain network topology (Bullmore and Sporns, 2012, Zhou et al., 2020) implies that long range connections between functional modules are more costly, and thus increased inter-modular connectivity should only be present in response to heightened computational demands. This theory is supported by studies which show that brain networks exhibit higher modularity and lower global efficiency at rest compared to when a task is being performed (Di et al., 2013, Kitzbichler et al., 2011), as well as during sleep compared to the awake state (Boly, 2012).

Past studies in TBI have shown that measures of functional segregation can be predictive of clinical and cognitive symptoms (Caeyenberghs et al., 2017, Kuceyeski et al., 2016). For example, Messé et al. (2013) demonstrated that modularity was inversely correlated with post-concussive symptoms in a cohort of individuals with mTBI at 6 months post injury, with higher modularity predicting more severe symptomatology. Similarly, Caeyenberghs et al. (2012) showed that individuals with severe TBI with higher connectivity degree and clustering coefficient displayed lower switching performance and more severe TBI symptoms. We did not find any significant correlations between symptom measures and graph metrics, possibly due to the late stage and mild nature of TBI in our cohort (1.88 total score out of a maximum of 12 on the TBI severity ranking).

We did find significant correlations between participant age and several of our graph metrics. There was a negative correlation between age and both global efficacy and clustering coefficient (these metrics decreased with older age). With modularity we saw the opposite trend, with higher baseline modularity being associated with aging. While we did not find a significant interaction between age and group (TBI vs. HC) for any network metric, these findings are in line with literature which has shown that network properties such as global efficiency and clustering coefficient are strongly correlated with aging and may predict the onset of dementia (Behfar, 2020, Sala-Llonch, 2014).

The present study does have some limitations that are worth discussing. The fact that this is the first investigation of its kind into chronic mTBI makes it difficult to compare our findings with other studies, limiting the context that can be used to facilitate a thorough interpretation of the data. It is important to note that the severity of mTBI in our cohort was on the mild end of the spectrum, and it is known that there are significant pathophysiological differences even within mild forms of TBI (Ettenhofer and Abeles, 2009, Konrad, 2011). Had our sample been more impaired, it is possible that we would have seen significant associations between clinical symptomatology and changes in network topology. Even so, few reports to date have directly correlated graph metrics with behavioral and clinical measures in TBI (Caeyenberghs et al., 2017). While there is some emerging evidence that certain network metrics, such as modularity, are associated with post concussive symptomatology (Caeyenberghs et al., 2017, Messé et al., 2013), the reported correlations are weak to moderate, and there are still other reports which found no relationships between graph metrics and clinical symptom scores (Kuceyeski et al., 2016, Hillary et al., 2014). Further work, perhaps using larger sample sizes and more advanced statistical analyses, are required to uncover what may be a complex, non-linear relationship between changes in network connectivity and behavior.

The heterogeneous nature of graph theory analysis methods is also an important source of variability in studies analyzing changes in network topology. The validity of any graph-based model of a complex system depends on the extent to which its nodes and edges represent true networks of the system under investigation. For some use cases, defining nodes and links is intuitive and straightforward. In social networks for example, nodes represent people and links are simply the connections between them (Lewis et al., 2008). In the case of brain networks however, defining nodes and links is much more complex (Fornito et al., 2013). Researchers can justify using a number of different parcellation schemes to define nodes, and methods for selecting significant connections can also vary widely (Sporns, 2018, Fornito et al., 2013). Small methodological differences can often yield divergent findings, making it difficult to compare results across studies. Compounding this issue is the fact that blood oxygen level-dependent (BOLD) signals are highly sensitive to noise artifacts, and that motion related noise can introduce a high degree of artificial correlation into graph-based metrics. For instance, in our study we found a significant main effect of the number of flagged volumes (a marker for image quality) on each of our outcome measures. While found no significant differences in the number of flagged volumes between groups or across time, it remains unclear the degree to which noise artifacts influenced the observed effects. Maintaining high data quality standards and using advanced preprocessing streams are important steps in mitigating the issue of nose contamination. Further sophistication of these methods would serve to improve the utility of network metrics as clinical and research tools. Despite its drawbacks however, graph theory analysis remains one of the most powerful tools with which to examine neural topology and characterize pathophysiology at the level of brain networks.

CRediT authorship contribution statement

Elias Boroda: Methodology, Formal analysis, Data curation, Writing - original draft, Writing - review & editing, Visualization. Michael Armstrong: Conceptualization, Methodology, Investigation, Supervision, Project administration. Casey S. Gilmore: Conceptualization, Methodology, Investigation, Supervision, Writing - review & editing. Carrie Gentz: Investigation, Supervision. Alicia Fenske: Investigation, Supervision. Mark Fiecas: Formal analysis. Tim Hendrickson: Formal analysis, Investigation, Writing - review & editing. Donovan Roediger: Formal analysis, Investigation. Bryon Mueller: Formal analysis, Investigation, Writing - review & editing. Randy Kardon: Conceptualization, Methodology, Investigation, Supervision, Project administration. Kelvin Lim: Conceptualization, Methodology, Investigation, Supervision, Project administration.