Age-dependent white matter disruptions after military traumatic brain injury: Multivariate analysis results from ENIGMA brain injury.

Mild Traumatic brain injury (mTBI) is a signature wound in military personnel, and repetitive mTBI has been linked to age-related neurogenerative disorders that affect white matter (WM) in the brain. However, findings of injury to specific WM tracts have been variable and inconsistent. This may be due to the heterogeneity of mechanisms, etiology, and comorbid disorders related to mTBI. Non-negative matrix factorization (NMF) is a data-driven approach that detects covarying patterns (components) within high-dimensional data. We applied NMF to diffusion imaging data from military Veterans with and without a self-reported TBI history. NMF identified 12 independent components derived from fractional anisotropy (FA) in a large dataset (n = 1,475) gathered through the ENIGMA (Enhancing Neuroimaging Genetics through Meta-Analysis) Military Brain Injury working group. Regressions were used to examine TBI- and mTBI-related associations in NMF-derived components while adjusting for age, sex, post-traumatic stress disorder, depression, and data acquisition site/scanner. We found significantly stronger age-dependent effects of lower FA in Veterans with TBI than Veterans without in four components (q < 0.05), which are spatially unconstrained by traditionally defined WM tracts. One component, occupying the most peripheral location, exhibited significantly stronger age-dependent differences in Veterans with mTBI. We found NMF to be powerful and effective in detecting covarying patterns of FA associated with mTBI by applying standard parametric regression modeling. Our results highlight patterns of WM alteration that are differentially affected by TBI and mTBI in younger compared to older military Veterans.

INTRODUCTION

Traumatic brain injury (TBI) is a significant health concern worldwide, with approximately 69 million diagnoses per year worldwide (Dewan et al., 2018). Military personnel are particularly at risk for TBI due to repetitive exposure to blasts or explosions (Hoge et al., 2008). Between 2000 and 2019, over 400,000 United States military personnel were diagnosed with a TBI, and approximately 84% of these diagnoses were mild TBI (mTBI) cases (Defense Medical Surveillance System et al., 2020). TBI diagnosis and severity are defined by duration of loss of consciousness (LOC), posttraumatic amnesia (PTA), and alterations of consciousness (AOCs) at the time of the injury (Kay et al., 1993). Diagnosis of mTBI limits LOC to less than 30 min, and PTA and AOC of less than 24 hr. Moderate and severe TBI diagnosis involve longer durations and may include abnormalities on conventional computed tomography (CT) imaging.

Over the past decade, there have been significant advances in neuroimaging techniques to assess TBI, particularly mTBI. Diffusion weighted imaging (DWI), which detects the movement of water molecules within brain tissue (Le Bihan, 2003; Mukherjee, Berman, Chung, Hess, & Henry, 2008), has been employed to indirectly evaluate the organization of white matter (WM) integrity. Fractional anisotropy (FA) measures the predominant direction of water molecule diffusion to infer the orientation of the long axes of axons and the fiber integrity. FA is not completely robust as it can provide misleading information related to crossing, diverging, or kissing fibers (Glenn et al., 2016). However, higher FA has been linked to increased directionality, suggesting healthy WM, whereas lower FA can implicate compromised WM.

Repeated exposures to TBI accelerate age‐related brain changes (Esopenko & Levine, 2015; Goldstein et al., 2012; McKee & Robinson, 2014). Specifically, reductions in FA, which are a normal part of healthy aging (Yap et al., 2013), are accelerated by repetitive blast exposures (Trotter, Robinson, Milberg, McGlinchey, & Salat, 2015). However, most prior reports on exposed military personnel have yielded inconsistent findings on the WM regions affected by TBI (Asken, DeKosky, Clugston, Jaffee, & Bauer, 2018). This lack of consensus, particularly pertaining to blast‐related mTBI, may stem from the tremendous spatial and individual heterogeneity due to the variability of biomechanical parameters related to blasts, including the directions and magnitudes of concussive forces, the presence of nearby rigid surfaces, the presence of protective gear, and other factors (Tate et al., 2021). Inconclusive findings in the literature may be partially explained by the presence of small statistical effects that are detectable only in large samples (Open Science Collaboration, 2015). Thus, the magnitudes and spatial distributions of injuries, coupled with the complexity of military‐related TBIs, imply that the accurate spatial mapping of WM disruption has remained elusive (Dennis et al., 2020).

Recent, large‐scale studies have used region of interest (ROI) approaches to study military TBI. However, head injuries, especially in a military context, exhibit heterogeneity in mechanism (as highlighted above), etiology (Davenport, 2016), and comorbidities such as post‐traumatic stress disorder (PTSD) and depression that make TBI‐related effects harder to interpret (Dennis et al., 2019; van Velzen et al., 2020). Thus, a data‐driven approach, unconstrained by the potential spatial biases associated with factors like anatomic parcellations, may empower the identification of WM damage patterns that are undetectable by ROI approaches due to these heterogeneous effects. By contrast, whole‐brain voxel‐wise methods impose stringent corrections for multiple comparisons that may lead to type II errors and that may obscure individual spatial heterogeneity across affected WM voxels.

To address these inferential challenges,

We accessed data shared by nine research cohorts in the Military Brain Injury working group from the Enhancing Neuro‐Imaging Genetics through Meta‐Analysis (ENIGMA) Consortium (Dennis et al., 2020).

We applied a multivariate, hypothesis‐free method, called non‐negative matrix factorization (NMF; Sotiras, Resnick, & Davatzikos, 2015; Xie, Ho, & Vemuri, 2011), to identify complex patterns of covariation (components) in data with high inter‐individual spatial heterogeneity. NMF uses a parts‐based approach to data representation (unlike principal component analysis or independent component analysis), used widely to identify image features.

We used regression modeling to test null hypotheses pertaining to associations between the mean FA of each component and the presence of military‐related TBI while adjusting for potentially confounding demographic and clinical variables. To understand cross‐sectional age‐related effects on the components, we explored the interaction between TBI diagnosis and age. Due to the high rates of PTSD and depression in military samples, we also controlled for the effects of these diagnoses.

We conducted the same analysis exploring mild TBI effects while adjusting for potentially confounding demographic and clinical variables by excluding sites with moderate or severe TBI participants.

Lastly, we conducted a confirmatory analysis with a more homologous sample by removing two older cohorts who recruited Vietnam‐era Veterans to ensure our analysis was not driven by these two older samples. We adjusted for potentially confounding demographic and clinical variables.

MATERIALS AND METHODS

Participants

Clinical, demographic, and imaging data of 1,475 participants from nine different cohorts included 725 Veterans who endorsed at least one military‐related TBI event and 750 Veterans without a TBI history (Table 1). Cohorts included United States Veterans and active‐duty service members, as well as military personnel from the Netherlands. Seven cohorts included personnel from the Iraq and Afghanistan military operations, and two studies included Vietnam‐era Veterans. All sites except the ADNI/DoD cohort were restricted to military personnel with mild TBI. Only the ADNI/DoD cohort of 201 participants included moderate and severe TBI. Due to heterogeneity in injury severity this cohort was removed from analyses using model_2 and model_3 described later in the manuscript. There were 63 participants (15 non‐TBI, 46 TBI) missing either PTSD or depression diagnosis information. Diagnostic criteria can be found in Table S1. We accounted for missingness of clinical data with a simple imputation method of replacing the missing variable with the mean value of that variable in participants from the same diagnostic category (i.e., missing PTSD for a TBI participant was replace with the mean PTSD value in the TBI group). T‐tests and chi‐square tests were used to determine group differences between demographic data (i.e., age, sex, PTSD diagnoses, and depression diagnosis) between TBI and non‐TBI groups. All participants provided written informed consent approved by local institutional review or ethics boards.

TABLE 1

Site demographics

Cohort	N	Scanner number	TBI (%)	Age (mean)	Females (%)	PTSD a (%)	Depression a (%)
ADNI DoD	201	20	115 (57.2%)	69.3	1 (0.5%)	88 (45.1%)	44 (22.6%)
Duke	298	3	160 (53.7%)	41.8	77 (25.8%)	64 (26.1%)	50 (20.4%)
INTRuST	85	11	54 (63.5%)	39.4	16 (18.8%)	34 (41.5%)	26 (31.7%)
iSCORE	118	1	41 (34.7%)	35.7	16 (13.6%)	54 (46.2%)	43 (36.8%)
MEDVAMC	49	2	35 (71.4%)	35.6	6 (12.2%)	44 (95.7%)	36 (78.3%)
VA Minneapolis	220	2	169 (76.8%)	32.6	12 (5.5%)	71 (39.7%)	79 (44.1%)
Stanford	35	2	45 (100%)	44.6	6 (17.1%)	NA b	16 (55.2%)
UMC Utrecht	94	1	10 (10.6%)	35.6	0 (0%)	46 (50.5%)	25 (27.5%)
VETSA	375	2	106 (28.2%)	61.8	0 (0%)	37 (9.9%)	42 (11.3%)
Total	1,475	44	725 (49.2%)	48.1	134 (9.1%)	438 (33.0%)	345 (26.0%)

Abbreviations: ADNI DoD, Alzheimer's Disease Neuroimaging Initiative‐Department of Defense; INTRuST, Injury & Traumatic Stress; iSCORE, Imaging Support for the Study of Cognitive Rehabilitation Effectiveness; MEDVAMC, Michael E. DeBakey Veterans Affairs Medical Center, UMC Utrecht, University Medical Center Utrecht; VETSA, Vietnam Era Twin Study of Aging.

Percentages were calculated excluding participants missing clinical diagnoses.

Stanford cohort did not collect PTSD diagnoses.

Preprocessing of diffusion imaging data

Each cohort's curators/owners processed their data locally using the ENIGMA DTI pipeline (http://enigma.ini.usc.edu/protocols/dti-protocols/#eDTI) for tract‐based spatial statistics (TBSS v1.2; Smith et al., 2006). Table S2 lists acquisition parameters for individual cohorts. Diffusion data were registered to the ENIGMA‐DTI FA mask and WM skeleton. In brief the ENIGMA‐DTI QC protocol involved three major steps: (a) examination of non‐normal distribution of mean FA values per region (tract) and removal of outliers (>2.69 × SD), (b) examination of diffusion weighted scans for each subject for scanner artifact and excessive motion, (c) inspection of FA maps after registration for misregistration of images. Visual inspection was conducted to remove data with artifacts or signal drop out. Visual inspection was done again by a central site to ensure proper registration. The WM skeleton specifies the center of the longitudinal axis along each WM tract. The mean FA skeleton was created for each subject by thinning the mean aligned FA values. The subject‐specific skeletons were then thresholded to include voxels with FA values exceeding >0.2 to exclude voxels in the cerebral spinal fluid and gray matter (Smith, Kindlmann, & Jbabdi, 2009). This removed voxels with both low mean values and with high inter‐subject variability within each site (Smith et al., 2006) and to ensure the skeletons to do not reach the outer edge of the cortex which can aid in alignment (Smith et al., 2009). Moreover, several studies examining the optimal thresholding of FA have reported results that are consistent with our selection of 0.2 (Kunimatsu et al., 2004; Taoka et al., 2009). The subject‐specific skeletons were then projected onto the ENIGMA skeleton by choosing the maximum FA values perpendicular to local skeleton structure within standard space. This approach aims to mitigate inter‐subject alignment errors (Smith et al., 2006). Individual subjects' skeleton FA data were shared by each participating site with the central site at Duke University to conduct the remaining steps of the analysis pipeline.

Nonnegative matrix factorization (NMF) analysis

NMF is a data‐driven approach that detects covarying patterns (components) within high‐dimensional data which has previously been used across a number of neuroscience applications (Cohen & Rothblum, 1993; Lee & Seung, 1999; Paatero, Tapper, Aalto, & Kulmala, 1991; Yang & Michailidis, 2016). More recently NMF has been applied to neuroimaging analyses (Sotiras et al., 2015; Xie et al., 2011), but to our knowledge has not yet to be applied to TBI analyses. The MATLAB code for calculating NMF components is available at https://github.com/asotiras/brainparts. We used NMF to find covarying patterns of FA across all 1,475 participants. First, the data were organized into an m × n matrix X where each row corresponded to the set of m voxels from the skeleton map, each column represented the set of n subjects, and each matrix element (i, j) contained the FA value of the ith voxel and jth subject. The data matrix X was given as input to NMF, which approximated it as a product of two non‐negative matrices W and H. The m × c matrix W represents m voxels (rows) and c components (columns), where c is specified by the user. Each matrix element (i, j) contains the loading of the jth component on the ith voxel, which denotes the relative contribution of each voxel to a given component. The above method is illustrated with a schematic in Figure 1. Components are estimated by positively weighting variables that consistently covary across participants to yield highly specific and reproducible patterns. The c × n matrix H corresponds to c components (rows), n subjects (columns), and each matrix element (i, j) contains the coefficient of the ith component and jth subject. These subject‐specific coefficients indicate the contribution of each component in reconstructing the original FA map. In general matrix factorization can be solved readily if the non‐negative constraint is removed. However, in an application such as the present, FA values that make up components cannot be negative, they must be positive. Factorization using a non‐negative constraint provides a solution that fits our context of white matter FA. Rather than including both positive and negative weights, NMF's nonnegative constraint allows for increased interpretability through a parts‐based representation of the WM (Sotiras et al., 2015).

FIGURE 1

Schematic of non‐negative matrix factorization (NMF) method. The data matrix X was given as input to NMF, which approximated it as a product of two non‐negative matrices W and H. The m × c matrix W represents m voxels (rows) and c components (columns), where c is specified by the user. Each matrix element (i, j) contains the loading of the jth component on the ith voxel, which denotes the relative contribution of each voxel to a given component

To determine the optimal number of components, we examined multiple NMF solutions ranging from two to 50, in increments of two, and assessed the reconstruction error for each solution as the Frobenius norm of the difference between the original data in X and the NMF approximation. We additionally evaluated metrics from a split‐sample reproducibility analysis, which was performed by splitting the data into two halves of similar age and sex distribution and examining how the reproducibility of the solution varies with the number of estimated components. The reproducibility was quantified by measuring the overlap between the independently estimated components for the two splits, after having matched them using the Hungarian algorithm (Kuhn, 1955). The overlap was measured using the inner product, which uses values in the range [0, 1], with higher values corresponding to higher overlap. Higher reproducibility was measured by the increasing overlap between individual components estimated by NMF in each split‐half of the data (Sotiras et al., 2015).

Statistical analysis

We sought to quantify the statistical association between the presence of TBI and the configuration of FA patterns captured by NMF‐derived components with linear mixed effects regression modeling using the fitlme function in MATLAB (v9.5). This function fits data to a linear mixed‐effects regression model, which can include both fixed and random effects. We first explored the main effects of TBI diagnosis of our sample while adjusting for age, sex, PTSD diagnosis and depression diagnosis as fixed effects and site/scanner as a random effect (model_1). We then examined the interaction between age and TBI diagnosis while adjusting for sex, PTSD diagnosis, depression diagnosis and site/scanner. Next, we excluded sites who recruited participants with moderate or severe TBI and explored the same interaction between age and mild TBI diagnosis while adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner (model_2). Lastly, as a confirmatory analysis to ensure our results were not driven by two older cohorts, we excluded sites who recruited Vietnam‐era Veterans to constrain our sample to military personal who served in Afghanistan and Iraq. We examined the same interaction between age and mTBI diagnosis with this younger sample (model_3).

We controlled for type I errors by correcting for the multiple comparisons associated with each component in model_1 and model_2 using the false discovery rate (FDR) method of Benjamini and Hochberg (Benjamini & Hochberg, 1995) implemented in MATLAB (v9.5). Corrected FDR significance values are reported as q‐values of p < .012 for model_1 and p < .004 for model_2.

RESULTS

The optimal number of components was derived from 1,475 participants who ranged in age from 18 to 85 (mean age = 48.1, SD = 15.7), with the majority (91%) self‐reported as male. For model_1, TBI and non‐TBI groups differed significantly by age (t[1,473] = −5.07, p < .001), sex (χ 2(1) = 4.24, p = .039), PTSD diagnosis (χ 2(1) = 45.9, p < .001), and depression diagnosis (χ 2(1) = 44.4, p < .001), with more females, older participants and less clinical diagnoses in the non‐TBI group. Model_2 was constrained to sites who recruited mild TBI (n = 1,278) and ranged in age from 18 to 71 (mean age = 44.77, SD = 14.14), with the majority (89%) self‐reported as male. The mTBI and non‐mTBI groups did not differ in sex (χ 2(1) = 3.49, p = .062), but significantly differed by age (t[1272] = −7.96, p < .001), PTSD diagnosis (χ 2(1) = 53.2, p < .001), and depression diagnosis (χ 2(1) = 48.6, p < .001), with older participants and less clinical diagnoses in the non‐mTBI group. Lastly, to ensure our results were not driven by our older cohorts of Vietnam‐era Veterans, model_3 included 899 military personnel and ranged in age 18–71 (mean age = 37.66, SD = 10.4), with the majority (85%) self‐reported as male. The younger mTBI and non‐mTBI groups did not differ in age (t[1272] = −1.07, p = .287), but differed by sex (χ 2(1) = 15.9, p < .001, PTSD diagnosis (χ 2(1) = 23.7, p < .001), and depression diagnosis (χ 2(1) = 29.7, p < .001), with more females and less clinical diagnoses in the non‐mTBI group. Above reported PTSD diagnosis and depression diagnosis statistics do not include imputed data.

Optimal number of components

The optimal number of components was determined to be 12 components. As displayed in Figure S1a, reproducibility overall decreased with an increase in number of components, yet a peak is demonstrated at 12 components. This inflection with a high inner mean product can be interpreted as relatively stable with a higher probability of reproducible networks for this dataset. Similarly, in Figure S1b, the reconstruction error decreased with an increase in network resolution, however there is a peak at 12 components suggesting a local performance maximum. While the reconstruction error and reproducibility are not consistent across the specified range of components, it does demonstrate stability in the 12 identified components. Fewer than the optimal number of components (i.e., 12) may omit necessary information modeled by NMF; however, adding more than 12 components may model uninterpretable information such as noise in the data due to motion. The optimal components (Figure S2) are composed of voxels that span portions of several anatomically defined white matter tracts (Table S3).

Visually, the 12 components display a high degree of anterior–posterior, radial, and left–right hemispheric symmetry. Components 2, 3, and 5 occupy the most central locations of white matter. Component 2 is primarily comprised of the bilateral anterior thalamic radiation and posterior limb of the internal capsule. Both components 3 and 5 encompass the corpus callosum, but component 3 also includes forceps major while component 5 involves forceps minor. More peripherally, components 10 and 11 are located. Likewise, component 10 spans forceps minor, but further extends to bilateral anterior corona radiate and inferior fronto‐occipital fasciculus. Similarly, component 11 includes bilateral inferior fronto‐occipital fasciculus as well as bilateral inferior longitudinal fasciculus. Next, component 6 primarily occupies bilateral anterior thalamic radiation, superior corona radiata, and posterior corona radiata. Component 8 encompasses bilateral superior longitudinal fasciculus, while component 7 spans forceps minor and bilateral inferior fronto‐occipital fasciculus. Component 1 and component 4 are comprised mostly of bilateral inferior longitudinal fasciculus. Component 12 is visually the most peripheral spanning forceps minor and inferior fronto‐occipital fasciculus. Lastly, component 9 occupies the most inferior brain locations including bilateral corticospinal tract and inferior cerebellar peduncle. The number of component voxels for each white matter tract is listed in Table S3.

Main effect and interaction effect of TBI

There was no significant main effect of TBI in model_1 (Table 2 and Table S4). However, we found a significant interaction (Table 3 and Table S5) after FDR correction between TBI and age after adjusting for sex, PTSD diagnosis, depression diagnosis, and site in four components including component 1 (q = 0.036), component 4 (q = 0.036), component 6 (q = 0.036), and component 12 (q = 0.036). FA values were negatively associated with age (as expected), but the association between age and lower FA was significantly stronger in the TBI than in the non‐TBI group (Figure 2).

TABLE 2

Model 1: effect of TBI, controlling for age, sex, PTSD diagnosis, depression diagnosis, and site/scanner

Component number	Voxel number	Beta estimate	Standard error	T statistic	Degree of freedom	Confidence interval	p value
1	30,217	0.00028	0.00147	0.190	1,469	−0.003, 0.003	.849
2	8,729	−0.00138	0.00189	−0.732	1,469	−0.005, 0.002	.464
3	4,117	0.00145	0.00259	0.557	1,469	−0.004, 0.007	.577
4	12,615	−0.00040	0.00168	−0.241	1,469	−0.004, 0.003	.810
5	4,084	−0.00299	0.00327	−0.915	1,469	−0.009, 0.003	.360
6	13,525	−0.00098	0.00144	−0.680	1,469	−0.004, 0.002	.497
7	10,740	−0.00010	0.00203	−0.050	1,469	−0.004, 0.004	.960
8	7,862	−0.00055	0.00183	−0.301	1,469	−0.004, 0.003	.764
9	5,764	0.00008	0.00182	0.043	1,469	−0.003, 0.004	.966
10	5,714	−0.00168	0.00196	−0.859	1,469	−0.006, 0.002	.391
11	5,719	−0.00195	0.00217	−0.899	1,469	−0.006, 0.002	.369
12	8,053	0.00058	0.00185	0.314	1,469	−0.003, 0.004	.754

TABLE 3

Model 1: interaction of TBI and age, adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner

Component number	Voxel number	Beta estimate	Standard error	T statistic	Degree of freedom	Confidence interval	p value	q value
1	30,217	−0.00023	0.00009	−2.518	1,468	−0.0005, −0.0001	.012**	0.036**
2	8,729	−0.00025	0.00012	−2.198	1,468	−0.0004, −0.00002	.028*	0.066*
3	4,117	−0.00013	0.00016	−0.838	1,468	−0.0004, 0.0002	.402	0.439
4	12,615	−0.00028	0.00010	−2.757	1,468	−0.0005, −0.0001	.006**	0.036**
5	4,084	−0.00027	0.00020	−1.356	1,468	−0.0007, 0.0001	.175	0.210
6	13,525	−0.00023	0.00009	−2.597	1,468	−0.0004, −0.0001	.009**	0.036**
7	10,740	−0.00026	0.00012	−2.108	1,468	−0.0005, −0.00002	.035*	0.066*
8	7,862	−0.00023	0.00011	−2.012	1,468	−0.0004, −0.00001	.044*	0.066*
9	5,764	−0.00022	0.00011	−2.015	1,468	−0.0004, −0.00001	.044*	0.066*
10	5,714	−0.00009	0.00012	−0.725	1,468	−0.0003, 0.0001	.469	0.469
11	5,719	−0.00025	0.00013	−1.916	1,468	−0.0005, 0.00001	.056	0.075
12	8,053	−0.00032	0.00011	−2.82	1,468	−0.0005, −0.0001	.005**	0.036**

p < .05.

Survived FDR correction (p < .012).

FIGURE 2

Model 1: interaction between age and TBI. Four components showed a significant interaction between TBI status and age in the total sample (n = 1,475). Component 1, component 4, component 6, and component 12 all displayed age‐dependent effects of lower FA that were stronger in TBI than non‐TBI groups

Main effect and interaction effect of mild TBI (mTBI)

We explored an analysis, excluding the cohort (n = 201) who recruited moderate and severe TBI participants, to address the age and mTBI interaction in the remaining 1,274 participants. There was no significant main effect of mTBI in model_2 (Table 4 and Table S6) We identified a significant interaction (Table 5 and Table S7) after FDR correction between age and mTBI while adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner in component 12 (q = 0.048). FA values were negatively associated with age (as expected), but the association between age and lower FA was significantly stronger in the mTBI than in the non‐mTBI group (Figure 3a).

TABLE 4

Model 2: effect of mild TBI, controlling for age, sex, PTSD diagnosis, depression diagnosis, and site/scanner

Component number	Voxel number	Beta estimate	Standard error	T statistic	Degree of freedom	Confidence interval	p value
1	30,217	0.0004	0.00165	0.235	1,268	−0.0028, 0.0036	0.814
2	8,729	0.0008	0.00189	0.411	1,268	−0.0029, 0.0045	0.681
3	4,117	0.0013	0.00236	0.562	1,268	−0.0033, 0.0060	0.574
4	12,615	0.0006	0.00175	0.319	1,268	−0.0029, 0.0040	0.750
5	4,084	0.0007	0.00283	0.254	1,268	−0.0048, 0.0063	0.799
6	13,525	0.0001	0.00158	0.048	1,268	−0.0030, 0.0032	0.962
7	10,740	0.0007	0.00214	0.306	1,268	−0.0035, 0.0049	0.760
8	7,862	−0.0002	0.00163	−0.147	1,268	−0.0034, 0.0030	0.883
9	5,764	0.0004	0.00198	0.207	1,268	−0.0035, 0.0043	0.836
10	5,714	−0.0006	0.00205	−0.277	1,268	−0.0046, 0.0034	0.782
11	5,719	−0.0005	0.00223	0.015	1,268	−0.0043, 0.0044	0.988
12	8,053	0.0009	0.00197	0.458	1,268	−0.0030, 0.0048	0.647

TABLE 5

Model 2: interaction of mild TBI and age, adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner

Component number	Voxel number	Beta estimate	Standard error	T statistic	Degree of freedom	Confidence interval	p value	q value
1	30,217	−0.0003	0.00011	−2.515	1,267	−0.0005, −0.0001	.012*	0.072*
2	8,729	−0.0002	0.00013	−1.198	1,267	−0.0004, 0.0001	.231	0.308
3	4,117	−0.0003	0.00016	−1.75	1,267	−0.0006, 0.00003	.080	0.137
4	12,615	−0.0002	0.00012	−2.057	1,267	−0.0005, −0.00001	.040*	0.134*
5	4,084	−0.0001	0.00019	−0.669	1,267	−0.0005, 0.0002	.504	0.550
6	13,525	−0.0002	0.00011	−1.834	1,267	−0.0004, 0.00001	.067	0.134
7	10,740	−0.0003	0.00014	−1.896	1,267	−0.0006, 0.00001	.058	0.134
8	7,862	−0.0002	0.00011	−1.62	1,267	−0.0004, 0.00004	.105	0.158
9	5,764	−0.0002	0.00013	−1.838	1,267	−0.0005, 0.00002	.066	0.134
10	5,714	−0.00007	0.00014	−0.493	1,267	−0.0003, 0.0002	.622	0.622
11	5,719	−0.0002	0.00015	−1.107	1,267	−0.0005, 0.0001	.268	0.322
12	8,053	−0.0004	0.00013	−2.87	1,267	−0.0006, −0.0001	.004**	0.048**

p < .05.

Survived FDR correction (p < .004).

FIGURE 3

Model 2: Interaction between age and mild TBI. Component 12 showed a significant interaction between mTBI status and age. The interaction in A was assessed in the sample (n = 1,274) of cohorts who only recruited mild TBI. The interaction in B further included a sample (n = 899) of cohorts who only recruited military personnel who served in the Iraq and Afghanistan military operations

Interaction effect of mTBI excluding older veterans

Finally, we conducted a confirmatory analysis, excluding the two cohorts with Veterans from the Vietnam‐era (n = 576), in 899 participants. We identified a significant interaction (Table 6 and Table S8) between age and mTBI while adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner in component 12 (p = .012). FA values were negatively associated with age (as expected), but the association between age and lower FA was significantly stronger in the mTBI than in the non‐mTBI group (Figure 3b). Due to the confirmatory nature of this analysis, an FDR correction was not imposed.

TABLE 6

Model 3: interaction of mild TBI and age, adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner after the removal of older Vietnam‐era Veterans

Component number	Voxel number	Beta estimate	Standard error	T statistic	Degree of freedom	Confidence interval	p value
1	30,217	−0.0003	0.00019	−1.574	892	−0.0007, 0.0001	.116
2	8,729	−0.0001	0.00021	−0.668	892	−0.0006, 0.0003	.504
3	4,117	−0.0003	0.00027	−0.940	892	−0.0008, 0.0003	.347
4	12,615	−0.0002	0.00019	−1.108	892	−0.0006, 0.0002	.268
5	4,084	−0.0001	0.00032	−0.443	892	−0.0008, 0.0005	.658
6	13,525	−0.0003	0.00018	−1.517	892	−0.0006, 0.0001	.130
7	10,740	−0.0003	0.00024	−1.157	892	−0.0008, 0.0002	.247
8	7,862	−0.0002	0.00018	−1.240	892	−0.0005, 0.0001	.215
9	5,764	−0.0003	0.00022	−1.534	892	−0.0008, 0.0001	.125
10	5,714	−0.0001	0.00022	−0.550	892	−0.0006, 0.0003	.583
11	5,719	−0.0003	0.00026	−1.220	892	−0.0008, 0.0002	.223
12	8,053	−0.0006	0.00023	−2.531	892	−0.0010, −0.0001	.012*

p < .05.

DISCUSSION

The present study applied data‐driven multivariate analyses to diffusion imaging and clinical data accessed through the ENIGMA Military Brain Injury working group (Dennis et al., 2020). We identified covarying patterns of fluctuation in voxel FA values, referred to as components, that were associated with TBI diagnosis. In four components, age‐dependent effects in lower FA were significantly stronger in the TBI group than in the non‐TBI group while adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner. When our dataset was constrained to mTBI, the same age‐dependent difference in FA remained in one component and was further identified when our data included a more homologous sample of Veternas from the Iraq and Afghanistan military operations. Our findings indicate that mapping FA patterns in diffusion data, defined by data‐driven methods, is an effective means to ascertain the impact of mTBI on WM organization in a manner that is unconstrained by predefined neuroanatomical tracts.

Military‐related mTBI involves unique mechanisms of injury including exposure mechanisms, comorbidities, and biological factors that may produce inter‐individual and spatial heterogeneity across white matter. Extensive evidence supports that blast TBI encompasses a unique pathophysiology (Goeller, Wardlaw, Treichler, O'Bruba, & Weiss, 2012; Salzar, Treichler, Wardlaw, Weiss, & Goeller, 2017), and repetitive exposure, which is commonly experienced in military conflicts, appears to have a dose–response relationship on white matter (Taber et al., 2015). After a head injury, which is often coincident with combat trauma, military personnel are at an increased risk of developing PTSD and depression (Stein et al., 2019) that negatively impact white matter (Davenport, Lim, & Sponheim, 2015; Dennis et al., 2019; Matthews et al., 2011). NMF can overcome this resulting heterogeneity by comparing subject‐specific coefficients derived for each component (Matrix H), where components are defined with a multivariate data‐driven approach, highlighting its power for future mTBI analyses. However, this method should be used with caution because component identification may be biased by systematic data quality issues. We still identified age‐dependent differences in our NMF‐derived components even when adjusting for highly comorbid disorders like PTSD and depression. White matter microstructural disruptions, in the form of lower FA, are associated with depression in Veterans with co‐occurring PTSD and TBI compared with Veterans with depression alone (Isaac et al., 2015; Matthews et al., 2011). Similarly, FA reductions are found in patients with both TBI and PTSD compared to TBI alone (Lepage et al., 2018). Controlling for the variance of these comorbid disorders is informative in understanding the complex factors that influence the relationship between TBI and white matter organization in military personnel.

Studies of military mTBI have yielded inconsistent findings on which WM tracts are affected and whether FA is increased or decreased (Asken et al., 2018). Published studies report lower FA after remote mTBI (Davenport, Lim, Armstrong, & Sponheim, 2012; Petrie et al., 2014), lack of significant mTBI effects (Hayes, Miller, Lafleche, Salat, & Verfaellie, 2015; Jorge et al., 2012; Sorg et al., 2016), and elevated FA in a sample similar to the present (Dennis et al., 2018). These inconsistencies may be due to the variability in mechanism and etiology surrounding mTBI. NMF is well suited to tackling inter‐individual spatial heterogeneity because our components are not neuroanatomically constrained by tract‐based ROI definitions, and at the same time, our method avoids the excessive burden of multiple comparison correction that plagues whole‐brain voxel‐wise methods (e.g., TBSS). Instead, NMF identifies patterns of FA‐fluctuation within the data unrelated to any previously defined ROIs. Each component with significant findings was composed of portions of multiple neuroanatomical white matter tracts. Prior data‐driven studies, which avoid these barriers, have demonstrated widespread white matter damage in military personnel and Veterans with mTBI (Davenport et al., 2012; Miller, Hayes, Lafleche, Salat, & Verfaellie, 2016; Morey et al., 2013; Taber et al., 2015). NMF constrains the coefficients to positive values and uses a parts‐based representation, which may be more interpretable than components derived from principal component analysis or independent component analysis. Prior findings applying NMF to healthy subjects identified NMF‐derived patterns of gray matter structural covariance that differed anatomically, but aligned closely with functionally defined brain networks (Sotiras et al., 2015). Future studies of mTBI should continue to consider hypothesis‐generating approaches such as NMF to further discern effects of mTBI as well as correlate data‐derived patterns with post‐injury functional outcome data.

Upon visual inspection, qualitative features that characterize the spatial distribution of the 12 components within the brain are noteworthy (Figure S2). The 12 components display a high degree of left–right hemispheric symmetry, a similar type of symmetry has been identified previously with NMF in patterns of gray matter (Sotiras et al., 2017). Component 12, which remained significant in all three analyses, shows anterior–posterior symmetry that affects WM almost exclusively in the frontal and occipital poles. In addition, component 12, along with component 4 and component 1, were the most peripherally located. While component 4 and component 1 did not remain significant after multiple comparison correction in the model considering only mTBI, these three peripheral components trended to exhibit stronger age‐dependent associations for lower FA in military personnel with mTBI. However, component 6 was significant when Veterans with moderate to severe TBI were included in the model. Contrastingly, component 6 encompasses more centralized white matter tracts, including hippocampal regions. Such a pattern is consistent with findings of damage to deep central structures typically limited to more severe injuries (Wilde et al., 2007).

Component 12 was the only component that remained significant across all three models. It primarily spanned the forceps major of the splenium, bilateral inferior fronto‐occipital fasciculus (IFOF), and bilateral inferior longitudinal fasciculus (ILF). One case study of a lesion in the forceps major described persistent deficits in the manipulation of visuo‐spatial information and specific types of navigation (Tamura et al., 2007), while global alexia was observed in other patients with injuries to this tract (Binder & Mohr, 1992). Studies using electrostimulation to the left IFOF resulted in deficits in semantic processing, specifically semantic paraphasia (Duffau et al., 2005; Duffau, Gatignol, Mandonnet, Capelle, & Taillandier, 2008; Moritz‐Gasser, Herbet, & Duffau, 2013) as well as non‐verbal semantic disruptions (Moritz‐Gasser et al., 2013). In addition, lower FA in the left IFOF has been correlated with worse semantic processing performance (Han et al., 2013). Stimulation to the right IFOF resulted in spatial cognition disruptions. Studies in stroke patients with damage to the right hemisphere spanning the right IFOF found left spatial neglect (Karnath, Rorden, & Ticini, 2009; Urbanski et al., 2011). Lastly, damage to the left ILF has been correlated with orthographic processing (Wang et al., 2020), while case studies of damage to the right ILF have visual hypoemotionality deficits (Fischer et al., 2016). Therefore, FA disruption of component 12 may be associated with several cognitive processes including visual attention, spatial cognition, semantic processing, and the integration of visual information and emotionality. However, the present study does not currently have neuropsychological test data to test these hypotheses.

Our results provide cross‐sectional evidence that age‐related decline in WM integrity is greater in individuals who experience mTBI. Relatedly, a trained machine learning algorithm that estimated brain age from the gray and white matter intensities of T1‐weighted MRI found the brains of TBI participants appeared 5–6 years older than the participants' chronological age (Cole et al., 2015). Using a whole‐brain voxel‐wise method in military Veterans exposed to blast forces, Trotter et al. (2015) similarly reported greater FA reductions with advanced age compared to blast unexposed Veterans. Even more concerning is the reported link between the chronic effects of TBI and neurodegenerative disorders, which exhibit progressive atrophy on structural MRI and FA decline of white matter based on DTI (Esopenko & Levine, 2015; McKee & Robinson, 2014). This evidence heightens concerns about TBI exposure and neurodegenerative disorders typically associated with advanced age, such as Alzheimer's disease, Parkinson's disease, or amyotrophic lateral sclerosis.

It is possible that our finding of stronger age‐dependent associations of lower FA in TBI could be explained by a single cohort of older Veterans in our sample with moderate/severe TBI (ADNI‐DoD), while the other cohorts in our sample were restricted to mTBI. This is worth noting particularly when coupled with evidence of a link between increased TBI severity and neurodegenerative disease later in life (Crane et al., 2016; Plassman & Grafman, 2015). To explore this assumption, we conducted an analysis after excluding this older moderate/severe TBI cohort. Our results yielded the same age‐dependent associations in one component even after a multiple comparison correction, which exhibited lower FA for military personnel with mTBI compared to military personnel without mTBI. This effectively demonstrates our initial significant finding in component 12 using the total sample was not driven by moderate and severe TBI in our older cohort. To ensure this finding was also not confounded by military personnel from a different era and theater experienced by our older Vietnam‐era Veterans, we explored a confirmatory analysis excluding a second older cohort (VETSA). The same age‐dependent association of lower FA in component 12 remained regardless of the restriction in age, sample size, and military era.

While it is unclear why the younger participants' FA is greater in the TBI group than the non‐TBI group, some explanations may be ventured. It is possible that young military personnel have differences in FA prior to an exposure to a blast or impact event, which may influence whether they report symptoms sufficient to meet criteria and be diagnoses with mTBI. In addition to this theory, military personnel with lower FA could experience more difficulty recalling events leading to the inability to meet criteria for mTBI during a retrospective interview at a younger age. More likely, younger military personnel may exhibit higher FA for a period of time after injury due to their recovery process. Acute inflammation could contribute to this higher FA (Kimura‐Ohba et al., 2016; Yang et al., 1999); however, this scenario is also unlikely due to the remoteness in time between the mTBI‐related event and data collection. Regardless, this might offer insight into the inconsistent results in the literature, where some studies report higher FA and some report lower FA in mTBI. However, these hypotheses can only be tested in a longitudinal design where diffusion imaging is acquired pre‐ and post‐exposure. Such a design could further incorporate multiple post‐exposure scans to confirm that age‐related FA changes are indeed accelerated by TBI, which can only be inferred by the present cross‐sectional design.

The data for this analysis were combined across nine cohorts including data acquisition site/scanner as a random effect in a mixed effects analysis. We evaluated leading methods to harmonize our diffusion imaging data across sites/scanners with particular emphasis on ComBat (Fortin et al., 2017; Hatton et al., 2020). ComBat is a method for harmonization of data from multiple sites or scanners using an empirical Bayesian method. Unfortunately, ComBat performance has only been tested for harmonizing data with as few as 20 samples per site, but not smaller. Several sites/scanners in our dataset consisted of fewer than 20 participants. Excluding these small sites would result in a dramatic attrition of our overall sample size. The ENIGMA Brain Injury consortium is also working on procedures to harmonize TBI‐specific clinical, symptom, and exposure parameters such as time since injury, mechanism, combat‐specific vs deployment‐specific injuries, blast‐related vs impact‐related injuries, post‐concussive symptoms, and cognitive functioning. Such procedures are expected to remove site/scanner‐associated variance from the mean‐FA within components for improved attribution to TBI‐associated variance. Nevertheless, we were able to identify patterns of disrupted white matter regardless of differences in scanning protocols and clinical assessments across sites by taking advantage of mixed‐effects modeling to adjust for site/scanner effects. Indeed, previous reports from ENIGMA suggest that the results from a mixed‐effects mega‐analysis are comparable to performing data harmonization with ComBat (Radua et al., 2020).

Combining neuroimaging datasets accessed through ENIGMA has allowed us to increase our sample size and using an approach like NMF has allowed us to reduce high‐dimensional neuroimage data, which provides increased power to detect smaller effect sizes. For instance, the weakest effect size among our significant findings was d = 0.17. Sample size estimation reveals that detecting this effect size with 95% power requires approximately 200 participants (Faul, Erdfelder, Buchner, & Lang, 2009; Faul, Erdfelder, Lang, & Buchner, 2007), which has been rare in the TBI literature. Further, there is a significant push within the scientific community to produce replicable results for brain‐wide associations that require large consortia‐type sample sizes (Marek et al., 2020). Performing large‐scale analyses across sites may also add to the goal of precision medicine (Thompson et al., 2020). Ongoing and future efforts to combine neuroimaging data from multiple cohorts will help to unravel the complexities and resolve inconsistencies in the TBI literature by increasing statistical power to detect differences that may be missed or irreproducible with smaller, underpowered samples.

LIMITATIONS

While our results are informative, the study has some limitations. First, it may be argued that our components could have been derived from our non‐TBI sample in order not to bias the selection of voxels that were subsequently compared to our TBI and mTBI groups. Rather, we chose to derive our components from the entire dataset to ensure significant power for the data‐driven analysis and to eliminate the possibility of confounding factors associated with only our non‐TBI sample. The sample was obtained from nine cohorts that acquired diffusion‐imaging data from 44 different sites and/or scanners, some with few participants per site. It may be beneficial to remove sites with limited sample size or harmonize the data prior to analysis as discussed above (Fortin et al., 2017; Radua et al., 2020). DTI analysis also have their own technique limitations such as incorrect estimations of fiber directions due to crossing, diverging, or kissing fibers. Future analyses can adopt more sophisticated methods for data collection such as high angular resolution diffusion imaging, diffusion spectrum imaging, or Q‐ball imaging (Glenn et al., 2016; Soares, Marques, Alves, & Sousa, 2013). In addition, a major challenge of research in military mTBI is that most studies depend on participant recollection and self‐report regarding exposures and symptoms often incurred during particularly chaotic times. This is generally deemed less reliable than contemporaneous accounts recorded by observers in the field. Furthermore, our study focused primarily on remote mTBI, so our results may not apply to the consequences of acute or subacute mTBI. Consortium efforts are underway to increase the sample size of our datasets and harmonize future variables such as time since injury, number of injuries, current symptomology, and cognitive outcomes that will be explored in future analyses.

CONCLUSION

We used NMF, an unsupervised, data‐driven method, to capture patterns of fluctuation across the white matter in the brain in a sample of military personnel accessed through the ENIGMA military brain injury working group. We identified significantly stronger age‐dependent effects of lower FA in the military personnel with a TBI compared to a military non‐TBI group. When our sample was constrained to only mTBI, we continued to see the same age‐dependent effect of lower FA for military personnel with mTBI even when controlling for potentially confounding clinical diagnoses. These findings highlight the power of combining data to perform large‐scale analyses across sites. Use of data‐driven methods may help to uncover the heterogeneous patterns of white matter damage resulting from military‐related mTBI and may uncover injury patterns in the brain unrelated to our current understanding of brain structure.

Data S1 Supplementary Information.

Click here for additional data file.

Table S1 Clinical diagnosis criteria for each cohort

Table S2. DTI acquisition parameters for each cohort

Table S3. Number of component voxels per each white matter tract

Table S4. Fixed effect, random effect estimate, and model fit statistics for model 1: effect of TBI, controlling for age, sex, PTSD diagnosis, depression diagnosis, and site/scanner

Table S5. Fixed effect, random effect estimate, and model fit statistics for model 1: interaction of TBI and age, adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner

Table S6. Fixed effect, random effect estimate, and model fit statistics for model 2: effect of mild TBI, controlling for age, sex, PTSD diagnosis, depression diagnosis, and site/scanner

Table S7. Fixed effect, random effect estimate, and model fit statistics for model 2: interaction of mild TBI and age, adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner

Table S8. Fixed effect, random effect estimate, and model fit statistics for model 3: interaction of mild TBI and age, adjusting for sex, PTSD diagnosis, depression diagnosis, and site/scanner after the removal of older Vietnam‐era Veterans

Click here for additional data file.

Figure S1 Optimal component determination. The optimal number of components was assessed from 2 to 50 in steps of two. Twelve components were chosen from the peak observed from the split sample reproducibility analysis and a low reconstruction error

Click here for additional data file.

Figure S2 Twelve independent NMF‐derived components. The 12 independent NMF‐derived components are visualized overlaying the ENIGMA DTI skeleton mask

Click here for additional data file.

Figure S3 Model 1 diagnostics: interaction between age and TBI quantile–quantile plots generated from the residuals for each of the 12 components in model 1 help visual the distribution of the data. Significant components for this model, including component 1, component 4, component 6, and component 12 show a normal distribution

Click here for additional data file.