Functional and Structural Neuroimaging Correlates of Repetitive Low-Level Blast Exposure in Career Breachers

Abstract

Combat military and civilian law enforcement personnel may be exposed to repetitive low-intensity blast events during training and operations. Persons who use explosives to gain entry (i.e., breach) into buildings are known as “breachers” or dynamic entry personnel. Breachers operate under the guidance of established safety protocols, but despite these precautions, breachers who are exposed to low-level blast throughout their careers frequently report performance deficits and symptoms to healthcare providers. Although little is known about the etiology linking blast exposure to clinical symptoms in humans, animal studies demonstrate network-level changes in brain function, alterations in brain morphology, vascular and inflammatory changes, hearing loss, and even alterations in gene expression after repeated blast exposure. To explore whether similar effects occur in humans, we collected a comprehensive data battery from 20 experienced breachers exposed to blast throughout their careers and 14 military and law enforcement controls. This battery included neuropsychological assessments, blood biomarkers, and magnetic resonance imaging measures, including cortical thickness, diffusion tensor imaging of white matter, functional connectivity, and perfusion. To better understand the relationship between repetitive low-level blast exposure and behavioral and imaging differences in humans, we analyzed the data using similarity-driven multi-view linear reconstruction (SiMLR). SiMLR is specifically designed for multiple modality statistical integration using dimensionality-reduction techniques for studies with high-dimensional, yet sparse, data (i.e., low number of subjects and many data per subject). We identify significant group effects in these data spanning brain structure, function, and blood biomarkers.

Introduction

Traumatic brain injury (TBI), particularly mild TBI, resulting from exposure to improvised explosive devices (IEDs) may be linked to long-term post-concussive sequelae and neuropathology. ^1,2 These observations have fueled intensive research efforts to understand the underlying acute and cumulative injury mechanisms of blast exposure. A significant component of these efforts has been to assess the effects of repeated low-intensity blast exposure in populations, such as military and civilian law enforcement “breachers” who use explosives to gain (breach) entry to buildings. One of the earliest studies of breachers evaluated Marine Corps personnel participating in a breacher training course at Quantico, Virginia. Course instructors with extensive previous exposure to low-intensity blast events and students with significantly lower blast exposures were included in the analysis. Functional neurological alterations, specifically those involving high-memory demand tasks, were observed in course instructors whereas no such decrements were observed in students participating in the course. ³

These findings raised questions concerning whether cumulative neurological alterations may occur in breachers over a career of repetitive low-level blast exposures. Subsequent studies of breacher populations have therefore focused on experienced (career) breachers with exposure to hundreds of blast events over many years. ^4,5 Breachers work within safe standoff distances guided by occupational standards ⁶ that are based upon risk for rupturing the unprotected tympanic membrane in humans ⁷ as well as overt pathology of the lungs and gas-filled organs in blast-exposed animals, as described over 50 years ago in seminal studies by Bowen and colleagues. ⁸ Although these current guidelines prevent injuries of acute clinical significance after blast exposure, they may not prevent subtle or cumulative neurological damage.

Both animal and human research examining the effects of blast have observed multi-scale neurological effects, including structural changes, ⁹ network-level dysfunction, ¹⁰ vascular damage and inflammation, ¹¹ and alterations in gene expression. ⁴ Given the heterogeneous nature of these effects, a variety of measures have been used to explore changes imputable to blast exposure and to capture information across the above-described domains. Traditional analyses for multi-modal studies involve performing separate statistical analyses per modality. For the imaging component of these studies, voxel-wise statistical testing ¹² is often used followed by a multiple comparisons correction step (for which various techniques ^13,14 have been proposed to handle the unique aspects of imaging data). Formal integration, however, is difficult given that such frameworks are not readily conducive to intermodality statistical inference, including those studies involving both imaging and non-imaging data. In addition, statistical power is a necessary consideration given that the expected effect sizes are potentially small and the number of data values far exceeds the number of subjects.

In order to address both of these issues in an exploratory, hypothesis-generating analysis, we use the statistical framework known as SiMLR (symmetric multi-view linear reconstruction ¹⁵ ), which performs data integration and dimensionality reduction before significance testing. SiMLR directly stems from earlier eigenanatomy work, ^{16 –19} which permits data-driven extraction of signals that explain population variability. Here, this framework is used across several modalities ¹⁹ to detect indicators suggestive of blast overpressure exposure, given that any signal(s) associated with overpressure exposure may be small, incomplete, and/or noisy. This is perhaps the first instance of using an integrative statistical design to detect effects of repeated low-level blast at the group level, spanning several neural systems and research modalities.

A variety of neurobiological measurements were used in this study to investigate blast exposure in career breachers. These measures are outlined in further detail below. The focus of the current report is the neuroimaging data collected as a secondary end point in the original study. Although the data integration framework uses all modalities (including fluid biomarker data and neuropsychological assessment) to determine neuroimaging findings, specific findings related to these non-imaging measures will be reported separately.

Methods

All procedures were reviewed and approved by institutional review boards at the National Institute of Neurological Disorders and Stroke, the Naval Medical Research Center, and Walter Reed Army Institute of Research. Informed consent was provided by all participants, and all procedures were performed during a 4-day evaluation at the National Institutes of Health Clinical Center. A total of 20 current or previous military or civilian law enforcement breachers were recruited into the study. Breachers must have had at least 4 years of experience in breaching and be actively involved (at least annually) in breacher training or operations.

Alternatively, former breachers could qualify if they have been exposed to at least 400 breaching blasts within a career. A total of 14 controls were recruited and exhibited similar characteristics as the experimental group in terms of age, sex, and operational experience. Controls must have had at least 4 years of military operational or law enforcement experience, be actively involved in military or civilian law enforcement training or operations, and could not have been exposed to more than 40 individual blasts over a career. A history of moderate-to-severe brain injury with loss of consciousness >5 min, diagnosis of central nervous system (CNS) disorder, and cardiac, respiratory, or other medical conditions affecting cerebral metabolism were all criteria for exclusion from the study.

Imaging

Magnetic resonance imaging (MRI) was performed using a 3 Tesla Siemens mMR Biograph (Siemens, Erlangen, Germany) with a 32-channel head and neck coil. Structural MRI sequences included high-resolution T₁-weighted magnetization-prepared rapid gradient echo (MPRAGE) imaging (repetition time [TR] = 1900 msec, echo time [TE] = 2.53 msec, flip angle = 9 degrees, and voxel size = 1  $\times$  1  $\times$  1 mm ${}^3$ ), T₂-weighted sampling perfection with application optimized contrasts using different flip angle evolution (SPACE) imaging (TR = 3200 msec, TE = 280 msec, and voxel size = 0.5  $\times$  0.5  $\times$  1 mm ${}^3$ ), fluid attenuation inversion recovery SPACE imaging (TR = 6000 msec, TE = 386 msec, and voxel size = 1  $\times$  1  $\times$  1 mm ${}^3$ ), and diffusion tensor imaging (DTI; TR 8900 msec, TE 98 msec, and voxel size = 2.2  $\times$  2.2  $\times$  3 mm ${}^3$ , 34 diffusion directions, with 55 slices).

In addition to structural MRI sequences, resting-state (rs-fMRI) functional MRI (fMRI) blood oxygen level–dependent (BOLD) data were obtained using a gradient-echo echo planar imaging sequence (EPI; TR = 2000 msec, TE = 30 msec, flip angle = 90°, spacing between slices 3.66, 3.0 mm slice thickness, and matrix size 64 $\times$ 64) in a 10-min acquisition, which resulted in 300 volumes. Brain perfusion was also determined utilizing first-pass gadolinium perfusion imaging (EPI pulse sequence, TR = 2800 msec, TE = 20 msec, flip angle = 40°, and voxel size = 2.5 $\times$ 2.5 $\times$ 4 mm ${}^3$ ). Forty sets of 30 images were acquired over a period of ∼1.5 min after administration of 9.7 mL of gadolinium contrast agent.

Demographics, clinical history, and neuropsychological assessment

Interviews were conducted with participants to collect demographic information and clinical history. Screening neurological history and examination were performed by a board-certified neurologist (E.M.W.). Psychometric testing was performed to evaluate a series of cognitive domains.

Testing included: Automated Neuropsychological Assessment Metrics (ANAM), Beck Depression Inventory–2nd edition (BDI-II), Booklet Category Test (BCT), Brief Symptom Inventory 18 (BSI-18), California Verbal Learning Test–2nd edition (CVLT-II), Combat Exposure Checklist (CEC), Delis-Kaplan Executive Function System (DKEFS) Sorting and Verbal subtests, Frontal Systems Behavior Scale (FrSBe), Immediate Post-Concussion Assessment and Cognitive Test (ImPACT), Interpersonal Reactivity Index (IRI), Iowa Gambling Task (IGT), Medical Symptom Validity Test (MSVT), NEO-Five Factor Inventory (NEO-FFI), Paced Auditory Serial Addition Test (PASAT), Revised Self-Monitoring Scale (RSMS), Rey Complex Figure Test (RCFT), Rivermead Postconcussive Questionnaire (RPQ), Satisfaction with Life Scale (SWLS), Trail Making Test (TMT), Wechsler Abbreviated Scale of Intelligence II (WASI-II), Wechsler Memory Scale III Logical Memory subtests (WMS-III LM I & II), and Wechsler Test of Adult Reading (WTAR).

Additionally, subjects were administered the Post-Traumatic Checklist–Military version (PCL-M), which demonstrated that neither breachers nor controls exceed the recommended post-traumatic stress disorder (PTSD) cutoff of 44.

Blood biomarkers

Non-fasting blood samples were collected using serum separator tubes between around 9:00 am and 12:00 pm, before interviews and the battery of tests. Samples were processed within 1 h of blood sample collection, using standard protocols, ²⁰ and then stored at −80 ${}^{\circ }$ C. Samples were stored until all samples had been collected, at which time batch assays were conducted.

Exosome isolation from serum

Exosomes were isolated from 0.8 mL of frozen serum. After thawing the serum, samples were centrifuged at 3000g for 15 min at 4 ${}^{\circ }$ C to remove debris and the supernatant was transferred into a clean tube. ExoQuick solution (System Biosciences Inc., Mountain View, CA) was added to the samples according to the manufacturer's instructions. Samples were incubated for 1 h at 4 ${}^{\circ }$ C and centrifuged at 1500g for 30 min. After the centrifugation, the supernatant was aspirated and the pellet was resuspended in 0.8 mL of Dulbecco's calcium- and magnesium-free salt solution (Sigma-Aldrich, St. Louis, MO). From each sample, 100 $\mu$ L of the resuspended pellet were stored at −80 ${}^{\circ }$ C and later used to measure total exosome protein content. The remainder of each sample was used to isolate and analyze neuronal exosomes.

To enrich for neuronal-derived exosomes, each suspension was incubated for 1 h at 4 ${}^{\circ }$ C with mouse/antihuman CD171 (L1CAM) biotinylated antibody (clone 5G3; eBioscience, San Diego, CA) in 50 μL of 3% bovine serum albumin (BSA; 1:3.33 dilution of blocker BSA in a 10% solution in phosphate-buffered saline; ThermoFisherScientific, Inc., Rockford, IL), for a final antibody concentration of 1:100. Next, 30 $\mu$ L of Streptavidin–agarose Ultralink resin (ThermoFisherScientific) in 50 $\mu$ L of 3% BSA was added to each sample tube and incubated overnight with continuous mixing at 4 ${}^{\circ }$ C. Then, samples were transferred to Pierce^TM centrifuge columns (0.8 mL; ThermoFisherScientific) and centrifuged at 6000g for 5 min. The content trapped in the column was then washed three times using 500 uL of 1 $\times$ Tris-buffered saline (TBS). For each wash, after adding the TBS, columns were centrifuged again at 6000g for 5 min. After the last wash, each column received 100 $\mu$ L of Pierce Gentle Ag/Ab Elution Buffer (ThermoFisherScientific) and was centrifuged at 10,000 rpm for 5 min. The neuronal-exosome–enriched sample was collected and stored at −80 ${}^{\circ }$ C.

Protein quantification

After having been thawed, samples received M-PER mammalian protein extraction reagent to lyse exosomes (ThermoFisherScientific), containing three times the suggested concentrations of protease and phosphatase inhibitors (complete™ ULTRA Tablets; MilliporeSigma, Burlington, MA). These mixtures were used to measure biomarker concentrations, using a site-specific Simoa HD-1 analyzer (Quanterix, Lexington, MA), an ultrasensitive paramagnetic bead-based enzyme-linked immunosorbent assay, according to the manufacturer's protocol. Coefficient of variation values were no higher than 30% for all analytes. Samples were randomized over plates, with laboratory scientists blinded to participant groups. For each sample, total and neuronal-enriched exosomal samples were run in the sample plate. All assays were run in duplicate.

Data processing

All imaging data described in the previous section were further processed using methods described. ^15,21 Specifically, T₁- and T₂-weighted MRI were used to generate cortical thickness maps; diffusion-weighted MRI was used to produce fractional anisotropy (FA) and radial diffusivity (RD) scalar images; rs-fMRI data were summarized using amplitude of low-frequency fluctuation (ALFF), fractional ALFF (fALFF), and network correlation spatial maps; and perfusion maps were generated from first-pass gadolinium perfusion imaging. Each derived scalar image was normalized to the subject's corresponding T₁-weighted image. The mapping from each T₁-weighted image to a cohort-specific template ²² was generated so that all data could be normalized to the same space. Image processing pipelines are described in further detail below.

Cortical thickness

The cortical thickness estimation pipeline ²³ found in the Advanced Normalization Tools (ANTs) package (https://github.com/ANTsX/ANTs) was used to generate regional cortical thickness maps from the T₁- and T₂-weighted MRIs. The processing of each T₁- and T₂-weighted MR image pair includes brain extraction, ²⁴ N4 bias correction, ²⁵ and six-tissue (cerebrospinal fluid, gray matter, white matter, deep gray matter, brainstem, and cerebellum) segmentation, ²⁶ the latter of which were used for determining a voxel-wise estimate of cortical thickness using the DiReCT (diffeomorphic registration based cortical thickness) approach. ²⁷

Diffusion-weighted imaging

The MRtrix suite of tools (https://github.com/MRtrix3) for processing diffusion-weighted imaging (DWI) was used for DTI reconstruction and generation of scalar images, such as FA. Pre-processing included denoising ²⁸ and whole-brain extraction. ²⁹ The pre-processed DWI were then converted to tensor images, ³⁰ FA, and RD maps. The average DWI was normalized to the corresponding T₁-weighted MRI using ANTs registration. ³¹

Resting-state functional magnetic resonance imaging

rs-fMRI data were processed using methods reported previously in the literature ¹⁵ and implemented in ANTsR (https://github.com/ANTsX/ANTsR). Each subject's BOLD time series was motion corrected by rigid normalization to an iteratively computed time-series spatial mean. Transformation parameters were kept for use as nuisance parameters and computation of frame-wise displacement. ³² Further nuisance variable calculation included contributions from non-cortical tissue with the CompCor framework. ³³ Transformation and CompCor nuisance variables were residualized from the BOLD signal, which was bandpass filtered (0.01–0.10 Hz) using wavelet filtering. ³⁴ Finally, fALFF ³⁵ and network correlation maps ³⁶ were calculated. Canonical network coordinates were provided by the Power coordinate system. ³⁶

Cerebral blood flow from first pass perfusion imaging

The voxel-wise mean of the control images was generated and used for motion correction of all time-point volumes. Cerebral blood flow (CBF) images were generated from the motion-corrected time series using a singular value decomposition deconvolution technique. ³⁷ The arterial input function voxel mask was generated automatically ³⁸ and manually inspected to ensure reasonableness of the generated results. Last, in our group-level analysis, we analyze the relative CBF by dividing each individual subject's CBF image by the global mean of CBF in the brain.

Template and template mapping

The symmetric group-wise normalization (SyGN) strategy was used for normalizing all imaging data to a common template ³⁹ for statistical analysis. The ANTs-based symmetric normalization (SyN) algorithm ³¹ was used to determine the correspondence between each subject's T₁-weighted image and the group template. ²³ Other modalities (e.g., mean BOLD) were aligned to the T₁-weighted image using ANTs linear registration tools. ⁴⁰ Transformation from each modality to the group template is performed by concatenating transforms such that only a single interpolation is performed. These mappings are illustrated in Figure 1.

FIG. 1.

Illustration of the template-based processing performed for each subject. Each modality is used to produce one or more derived images (indicated by yellow arrows), which are then transformed to the group template for statistical analysis by each subject's T₁-weighted image. “↔” denotes a linear mapping whereas “” denotes a diffeomorphic deformable mapping. DWI, diffusion-weighted imaging; FA, fractional anisotropy; fALFF, fractional amplitude of low-frequency fluctuations; rs-fMRI, resting-state functional magnetic resonance imaging; T1-w, T₁-weighted; T2-w, T₂-weighted. CBF, cerebral blood flow.

Statistical analysis

Exposure to blast overpressure can lead to a cascade of effects that vary across scale and time. Alterations in brain structure and function, quantified through various imaging modalities, can cause metabolic, structural, and functional network-level dysfunction, thereby disrupting cognition—the latter observed through self-report, clinical visits, and/or neuropsychological assessment. Concomitant vascular and inflammatory effects can be assessed through the use of fluid biomarkers. Each signal in the cascade of effects may be small, incomplete, and/or noisy, but potentially mutually overlapping in contributing to our understanding of blast effects.

SiMLR ⁴¹ is a statistical technique that derives from earlier work, specifically prior-based eigenanatomy ^17,42 and sparse canonical correlation analysis in neuroimaging (SCCAN). ^18,43 These techniques invert the typical analysis paradigm which performs statistical testing before voxel-wise pruning and clustering to correct for multiple comparisons as well as enforce spatial smoothness in the context of image data. ¹² In contrast, eigenanatomy-based techniques cluster data and reduce dimensionality as a pre-processing step before hypothesis testing. This order of operations facilitates data inspection for outliers; visualization of correlations across the extent of the brain before tests for specific outcomes or group differences is performed.

Last, after the completion of data cleaning, one performs statistical testing on these data-derived clusters, similar to standard principal component regression, ⁴⁴ but for both imaging and tabular data. This is also similar in spirit to other popular frameworks such as region of interest–based analysis ⁴⁵ or tract-based spatial statistics ⁴⁶ with dimensionality reduction determined by specifying multiple anatomical regions, often on the order of 10–100 regions. However, this anatomical simplification is not purely data derived by some optimality constraint nor do they permit formal multi-modal integration. SiMLR, in contrast, strikes a balance between defining a priori regions and allowing data-driven exploratory studies of signals that are distributed across networks of variables in imaging and/or tabular formats (Fig. 2). Note that all these techniques are available as open source through the ANTs software ecosystem (https://github.com/ANTsX/ANTsR).

FIG. 2.

The various effects observed in traumatic brain injury overlap, but the signals used for assessing the various components are often noisy, incomplete, and/or small. SiMLR leverages the covariation between imperfect data sources (along with dimensionality reduction techniques) to infer the underlying signal(s) of interest in a focused, statistically powerful way. PCA, principal component analysis; SiMLR, similarity-driven multi-view linear reconstruction.

Whereas earlier SCCAN work was limited to examination of two modalities, SiMLR permits an arbitrary number of modalities and customized regularized models as input. Output consists of low-dimensional embeddings for each modality that best predict the other $N-1$ input modalities using corresponding bases derived from matrix decomposition techniques, such as singular value decomposition or independent components analysis. For the ith modality, X_i , this solution consists of a pair of matrices $\left(U_i,V_i\right)$ where the columns of U_i describe the low-dimensional bases spanning the space of X_i , where the number of bases is a user-tuned parameter. The columns of V are the regularized, sparse representation of each modality component. These solutions are iteratively determined by, first, reconstructing each modality from the basis functions of the other modalities, that is, $U_{j\ne i}$ , and, second, updating each modality's basis functions (Fig. 3). ⁴⁷

FIG. 3.

Illustration of inferential design using SiMLR. Multiple modalities are projected into a regularized, low-dimensional space for linear regression analysis (similar to PCA regression). In the imaging context, for a specific modality, each of these low-dimensional projections comprise a similarly varying (across subjects) spatial region in the space of the template. DMN, default mode network; FA, fractional anisotropy; fALFF, fractional amplitude of low-frequency fluctuations; PCA, principal component analysis; RD, radial diffusivity; SiMLR, similarity-driven multi-view linear reconstruction.

For each imaging modality, this decomposition results in spatially coherent regions grouped by variation similarity across subjects in the space of the template. Table 1 summarizes the data utilized in SiMLR. Subsequent statistical analysis permits determination of statistically significant regions and visualization of these physical locations within the brain template. SiMLR processing is unsupervised and produces a low-dimensional space that jointly spans the neuroimaging data and other included modalities. Subsequent to the SiMLR-based decomposition, multiple linear regressions were performed in the low-dimensional space for each modality basis. Specifically, we perform regressions of the form:

Table 1.

Summary of the Various Data Used in the SiMLR Statistical Framework

Demographics • Age • Education • Number of breaches (last year) • Number of breaches (career)	Neuroimaging • Structural ○ Cortical thickness (by T1-/T2-w MRI) ○ FA (by DWI) ○ Radial diffusivity (by DWI) • Functional ○ fALFF (by rs-fMRI) ○ default mode network (by rs-fMRI) ○ β (perfusion, by gadolinium-enhanced contrast)
Psychometric batteries • Wechsler Test of Adult Reading • Wechsler Memory Scale • Wechsler Abbreviated Scale of Intelligence • Rey-Osterrieth Complex Figure Test • Immediate Post Concussion Assessment and Cognitive Testing • Automated Neuropsychological Assessment Metrics (ANAM4) • Booklet Category Test • Delis Kaplan Executive Functioning- Sorting Test	Fluid biomarkers (exosome) • Inflammation ○ Interleukin-6 ○ Interleukin-10 ○ Tumor necrosis factor-α • Neural injury ○ Neurofilament light ○ Tau ○ Aβ42

Although not all data collected were used, a selection encompassing the various aspects of effects observed in traumatic brain injury were selected.

SiMLR, similarity-driven multi-view linear reconstruction; T1-w, T₁-weighted; T2-w, T₂-weighted; MRI, magnetic resonance imaging; FA, fractional anisotropy; DWI, diffusion-weighted imaging; fALFF, fractional amplitude of low-frequency fluctuation; rs-fMRI, resting-state functional magnetic resonance imaging; Aβ42, beta-amyloid 42.

Table 2.

Demographic Summary of the Experienced Breacher and Control Cohorts

	Breacher (N = 20)	Control (N = 14)
Age, years	39.65 (26–54)	38.85 (27–53)
Ethnicity	• Caucasian = 17 • American Indian = 1 • Asian = 1 • Other = 1	• Caucasian = 12 • African American = 1 • Asian = 1
Handedness	Right = 18; left = 2	Right = 12; left = 2
Service	• Army = 13 • Navy = 1 • Law enforcement = 6	• Army = 10 • Law enforcement = 4
Duration of service	16.8 (8–27)	13.8 (5–24)
Post-Traumatic Checklist–Military (PCL-M)	20.7 ± 4.7 (17–30) with 0 subjects >44 (PTSD cutoff)	25.6 ± 7.1 (17–40) with 0 subjects >44 (PTSD cutoff)
Concussion	Yes = 13; no = 7	Yes = 6; no = 8
Total blast exposures	• 4628 (100–34,800)	3 (0–35)
Breaches in past year	• >400 = 3 • 200–399 = 5 • 100–199 = 3 • 40–99 = 1 • 10–39 = 2 • 1–9 = 4	• None = 14

PTSD, post-traumatic stress disorder.

$$<italic>U_{i,l}\sim Age+Education+BrainVolume+Group<∕italic>$$

where i denotes the modality and l denotes the $l^{th}$ basis vector. The brain volume for each subject was calculated directly from the brain extraction step of the cortical thickness processing. Based on the limited number of subjects and previous experience, we specified two basis functions per modality.

Results

Twenty career breachers were recruited into the study and compared to 14 matched controls. The original study was designed for 20 controls, and it was estimated that 16 persons per group were required to detect a difference between the breacher and control groups on the primary outcome measure, memory performance on the Automated Neuropsychological Assessment Battery. ³ Fourteen military and 6 law enforcement persons were included in the breacher group, whereas the control group comprised 10 military and 4 law enforcement personnel. No differences were observed between groups with respect to age, ethnicity, handedness, service, duration of service, or report of previous concussion. Career breachers reported an average of 4628 (100–34,800) breaching blast exposures over their careers, whereas controls reported an average of three (0–35) exposures. Eighteen of 20 career breachers reported exposure to breaching blast within the past year, whereas none of the controls reported exposures in the previous year.

Using the SiMLR methodology, statistically significant differences were observed between breachers and controls on structural and functional neuroimaging. In Figures 4–6, we visualize the network-like regions derived from SiMLR for each of the imaging modalities in the space of the template. Color maps are chosen to represent each individual modality's weighted average score, not to imply greater or lesser importance of any individual voxel within the colored region. The exception, here, is the default mode network (DMN) where we explicitly show that the weighted feature vector does not overlap with the default network and, as such, we use two color maps. Specifically, Figure 4 showcases cortical thickness and perfusion results. Figure 5 illustrates DWI results using FA and RD maps. Finally, functional activity by fALFF and default mode internetwork correlation maps are provided.

FIG. 4.

Group differences in cortical thickness (top) and perfusion (bottom) imaging. On the left are box plots showing the subject distribution based on the SiMLR-based imaging-specific decompositions. Along the x-axis are the control versus breacher groups with the raw p value in the upper left of the plots. Cortical thickness projections showed a general cortical thickening (driven largely by occipital lobe and default mode network regions). No group differences in relative perfusion survived multiple comparisons correction. SiMLR, similarity-driven multi-view linear reconstruction.

FIG. 5.

White matter group differences illustrated with FA and RD projections derived from DWI. On the left are box plots showing the subject distribution based on the SiMLR-based imaging-specific decompositions. Along the x-axis are the controls and breacher groups with the raw p values in the upper left of the plots. Whereas fractional anisotropy is decreased in the breacher group (vs. controls), radial diffusivity is increased. Both measurements are driven primarily by frontal lobe white matter. DWI, diffusion-weighted imaging; FA, fractional anisotropy; RD, radial diffusivity; SiMLR, similarity-driven multi-view linear reconstruction.

FIG. 6.

Functional connectivity group differences illustrated with DMN correlation maps and fALFF projections derived from resting-state fMRI. On the left are box plots showing the subject distribution based on the SiMLR-based imaging-specific decompositions. Along the x-axis are the controls and breacher groups with raw p values in the upper left of the plots. Connectivity between default mode network and other networks is increased in breachers whereas resting-state activity, largely within the default mode regions, is higher in controls. DMN, default mode network; fALFF, fractional amplitude of low-frequency fluctuations; fMRI, functional magnetic resonance imaging; SiMLR, similarity-driven multi-view linear reconstruction.

In addition, to further summarize the overall pattern of differences across modalities, we provide a spider plot of the SiMLR projections between control and breacher groups, in Figure 7, across all measurements that survive family-wise error rate correction. We used Bonferroni adjustment by multiplying all raw p values by 16. This significance cutoff corresponds to a p value of 0.003125. All results that survive this correction demonstrate p less than or equal to that value. Further, we note that p values are not surrogate measurements for the strength of an effect, and we do not propose to quantify relative importance of predictors within the set that survive correction.

FIG. 7.

Summary statistical findings, in terms of mean modality-specific projection values, over all modalities used in the study, demonstrating statistically significant differences between breachers and controls. Eight modalities × two bases per modality resulted in 16 statistical tests. All axes of the spider plot represent significant effects where the limits of each axis represent the minimum and maximum of the modality-specific projections. Blood biomarker results, though important to this study, will be detailed in a separate publication where greater specificity about the role of each exosome measurement will be reported. DMN, default mode network.

Shown in Figure 4 are the group structural imaging differences in both cortical thickness (top) and gadolinium-based perfusion (bottom). Specifically, the sparse representation reveals a pattern of relatively increased cortical thickness in experienced breachers versus controls. The primary neuroanatomy involved in the cortical component includes the DMNs as well as visuospatial regions in the occipital lobe. The perfusion component, on the other hand, does not show a group-level difference that survives Bonferroni correction. However, breachers exhibit a trend toward reduction in relative perfusion, focusing in the superior frontal gyrus, caudate nuclei, and posterior putamen. Scatter plots, along with the mean and 95% confidence intervals, for each group are shown for the corresponding projection data for each modality. Only thickness results are statistically significant at the Bonferroni level ( $p<0.0016$ ; $q<0.0257$ ).

RD and FA results are found in Figure 5. Group differences in both measurement types, derived from the diffusion tensor at each voxel, are largely concordant. Both components include the corpus callosum (CC) with stronger weights (signified by brighter voxels) in the anterior extent of the CC (the genu) in the RD component. Portions of the claustrum and splenium contribute to these effects in the FA component. Effects are significant at the $p<0.0026$ level ( $q<0.041$ ).

Functional differences are illustrated in Figure 6 for both the DMN and fALFF measurements. Internetwork connectivity from the DMN to normally anticorrelated task-control regions is higher in breachers. Further, higher resting activity is observed in controls with the fALFF modality, and this effect overlaps substantially with default mode regions. Both sets of differences are statistically significant ( $p<2e-4$ ). There was no evidence of motion differences between breachers and controls ( $p<0.5$ ). The modalities/measurements that survive Bonferroni correction are summarized in the spider plot shown in Figure 7. Components of structural (cortical thickness, RD, and FA) and functional (default mode network and fALFF) imaging all demonstrate significant effects after correction, as do blood biomarkers. Neuropsychometric scores did not survive Bonferroni or false discovery correction in this study (p > 0.05, q > 0.05). Neuropsychological assessment and blood measurement results will be detailed in separate publications.

Last, blood biomarker results, though important to this study and showing significant relationship with breacher group status, will be detailed in a separate publication. In this forthcoming study, greater specificity than what is available from SiMLR about the role of each exosome measurement will be reported.

Discussion

This study was conducted to address the need for an improved understanding of potential cumulative neurological sequelae after repetitive low-level blast exposure. It is the first study of experienced breachers with a history of repeated blast exposure to include a wide spectrum of behavioral and imaging modalities. This work aligns with guidance established in fiscal years 2018 and 2020 United States National Defense Authorization Acts (NDAA; Public Law Nos. 115-91 and 116-92) to better understand the health and performance effects of blast exposure on members of the Armed Forces during combat and training.

Many previous studies on the neurological effects of blast exposure have focused upon clinically diagnosed TBI after a significant blast event, such as exposure to an IED. ^{48 –51} Of note, these types of exposures are often accompanied by other, better studied, traumatic neurological insults such as acceleration-deceleration, impact, or shrapnel-related injury, often resulting in polytrauma rather than isolated blast injury to the CNS. Whether this cumulative exposure results in diagnosable injury remains unclear. The current study was conducted to better understand whether career breachers manifest neurological findings that could represent cumulative injury and to identify potential areas for future hypothesis-based studies. Breachers are a unique population for further study, given that their exposure to blast often occurs in the absence of other traumatic neurological insults, thus allowing for a more direct study of the neurological sequelae of primary blast overpressure exposure.

This study is also the first to analyze several modalities in an integrative statistical framework for detecting the effects of chronic blast exposure. Multiple imaging modalities were used to capture a diverse set of measurements of the structure and function of each subject within the cohort. These imaging data, coupled with neuropsychological assessments and exosomal data, permitted a focused, data-driven approach to identifying variables which distinguish controls from breachers, despite the limited number of subjects, the large data dimensionality, and related concerns regarding statistical power.

We found a general increase in cortical thickness throughout the cortex, especially within occipital lobes, in career breachers as compared to controls. Most earlier studies exploring cortical measures in TBI, including those exploring blast-related TBI, have demonstrated cortical thinning in injured persons compared to controls, ^{9,52 –55} whereas others have exhibited increases in thickness. ^56,57 However, in contrast to the current study, these earlier reports involved persons with clinically diagnosed injury after blast exposures, many of which were accompanied by conventional neurological insults as well. Therefore, the significance of the differences in this study is uncertain. The finding of cortical thickening in career breachers accompanied by concomitant increases in RD and decreases in FA-DTI measures could reflect reduced cortical myelination or, potentially, less organized intracortical connections.

Alternatively, these findings might be related to astroglial scarring at boundaries between gray matter and white matter, as recently described in a post-mortem case series of human brain after blast exposure. ⁵⁸ Alterations in cortical myelination, intracortical connections, or glial scarring at the gray-white interface could all result in an apparent increase of the cortex metric using pipelines that are designed to distinguish between gray and white matter in the calculation of this measure.

Given its sensitivity to structural changes and reflection of physiological alterations in the brain, perfusion studies have been increasingly used in studying TBI. ^{59 –62} Vascular measures are of particular interest in the current effort given mounting pre-clinical evidence of altered blood–brain barrier integrity, vascular reactivity, and endothelial ultrastructure after blast overpressure exposure. ^{11,63 –65} Career breachers demonstrated an overall decrease in relative CBF (perfusion) in comparison to controls, particularly in the caudate nucleus, somewhat consistent with previous studies. ⁶⁶ From the current data, it is unclear whether these findings reflect fixed perfusion defects attributable to tissue-level structural alterations, changes in vascular reactivity, and/or alterations of metabolic demand within regions demonstrating reductions in CBF.

Diffusion-based imaging acquired in multiple directions with the production of scalar images from DTI-based measures such as RD or FA can provide insight into the integrity of white matter tracts throughout the brain. In this study, we demonstrated a simultaneous decrease in FA and increase in RD in breachers versus controls within the genu and posterior portions of the CC, along with corticospinal tract and claustrum. These findings are consistent with reductions in white matter integrity within these regions. No significant differences in motion or motion parameters existed between groups in the DWI data. Earlier work exploring DTI-based measures in TBI is extensive. ^{67 –71} Whereas FA and RD may demonstrate increased, decreased, or bidirectional values in the acute to subacute phase post-TBI, ⁷² chronic white matter injury largely presents with decreased FA and increased RD, ^68,69,73 consistent with our findings.

Specific to studies involving white matter integrity after blast exposure, several groups have demonstrated alterations using various DTI-based measures across different time scales. ^69,74,75 However, in contrast to the current effort, these previous studies have focused upon clinically diagnosed blast-related TBI, which may involve conventional insults (e.g., impact, acceleration-deceleration) in addition to primary blast overpressure. Further, many of these studies are strictly looking at one or more DTI-based measurements without the benefit of leveraging other imaging and non-imaging data sources in a joint analytical framework.

Finally, functional differences were found in the connectivity between default mode (DMN) and other networks, as well as measurements of resting-state activity (fALFF). Specifically, increased connectivity between DMN and other networks and reductions in DMN activation were observed within breachers compared to controls. Functional connectivity analysis is particularly sensitive to the effects of motion, which produce spurious correlations. Critically, we showed that no significant differences in motion or motion parameters existed between groups. In TBI, increases in connectivity between DMN and other networks have been previously reported. ^76,77 Additionally, reductions in DMN activation after TBI have been described as well ^78,79 and are consistent with observations in the current study.

Last, blood measurements and neuropsychological assessments (at trend level) suggest that these changes detected on imaging may be associated with other differences between groups. Exosomal measurements, as used here, suggest increases in neuroinflammation (by interleukin [IL]6, IL10) and neural injury (by tau, ab42, and neurofilament light). Whereas these modalities are not the primary focus of this work, they are consistent with the theory that exposure to breaching-related blasts leads to system-wide effects in the brain.

Although the imaging findings are statistically robust and consistent with some potential mechanisms of injury, their attribution to blast exposure has to be considered provisional at this point. Imaging was included in the original study as an exploratory outcome, and the sample size was not based on an estimate of the statistical power required to detect a difference between the two groups on any single comparison. Moreover, an additional limitation of this study is that breachers have slightly higher PCLM scores than controls. Although neither breacher nor control subjects exhibit PCLM scores that meet a standard PTSD cutoff of 44, ⁸⁰ future work, with greater variability in the presence of PTSD or more optimally balanced groups, may be able to separate effects of blast exposure from PTSD. However, power in this cohort is too limited to explore this question.

Relatedly, although these subjects were richly quantified, statistical power prevents us from fully testing every measurement in every modality, which would result in, potentially, hundreds to millions of tests, depending on the image resolution used or the collection of brain regions used. Bonferroni correction over this number of tests would be prohibitive. As such, we used SiMLR to compute a reduced dimensionality space that permitted a relatively small number of inferential tests to be performed, here two for each modality (16 total with Bonferroni cutoff of p < 0.003125). SiMLR allows us to retain a degree of interpretability, relative to a method such as singular value decomposition, in that the derived components are spatially constrained to brain tissue, are non-negative, and are sparse, that is, select for only a portion of the brain per component in a data-driven manner.

Although SiMLR provides a sparse representation for each of the derived variables—which is suggestive of a univariate statistical map—Figures 4 through 6 must be interpreted with caution in that they are fundamentally multi-variate: The outcome variables derived from these maps involve each part of the brain shown in the figures, and no individual voxel can be treated as independently more important than the others. Although this is appropriate for network-like interpretations, we lose a degree of spatial specificity in that we cannot isolate an effect to a single structure. As always, there is a balance between interpretability, specificity, and statistical power when performing inferential studies, and the balance struck in this work is not the only possible valid approach.

In summary, we found statistically significant differences on a number of imaging measures between a group with extensive occupational blast exposure and unexposed controls. Focused, hypothesis-driven studies will be required to establish the reproducibility of these changes and attribute them with confidence to the effect of blast. Additional work is needed to look more focally at cortical and white matter regions, determine a potential set of causal pathways that underlie reported performance deficits in career breachers, and study in more detail whether mechanisms that link the signals observed from exosomal data to multiple MRI modalities. Further, these findings should be extended into other populations, such as artillery and explosive ordinance disposal service members. Ultimately, additional research cohorts will aid in refining safety protocols for military and civilian personnel routinely exposed to blast. Additional work may also help to determine factors that could confer resilience or vulnerability to low-level blast exposure at the individual level.