Changes in volumetric and metabolic parameters relate to differences in exposure to sub-concussive head impacts

Participants

This study was approved by the Queen’s University Health Sciences Research Ethics Board (Kingston, ON, Canada), and informed consent was obtained from each participant. The study number was SURG-399-17 (file number: 6021022). The study was conducted in accordance with the Declaration of Helsinki. A total of 44 Canadian male collegiate football student-athletes, varying in positions, were recruited to participate in this study. All subjects were scanned before the season (‘PRE’ group), within two months prior to participation in full contact practices (Figure 1). None of the athletes enrolled in this study were recovering or suffering from post-concussion symptoms at the time of the scan, or in the year prior. OPEN IN VIEWER

Neuroimaging was repeated in-season following the 14-day training camp period and the first two regular season games (post-training camp; ‘PTC ’ group), and post-season, one month following the last competitive game (‘POST ’ group; Figure 1). Of the 44 players recruited from baseline testing, nine were removed from the PTC cohort prior to scanning (three concussions and six season-ending injuries), for a total of 35 follow-up subjects. Of these 35 athletes, two sustained a concussion during the second half of the season, and two suffered a season-ending injury. Thus, a total of 31 players (Figure 1) from the original cohort returned for the final time point (POST), approximately 30 days after their last participation in contact activities.

At each time point, subjects were also asked to complete the symptom evaluation section of the Sport Concussion Assessment Tool (3rd edition; SCAT-3), ²³ in order to monitor changes in self-rated symptoms throughout the season.

Monitoring of head impact exposure per session

Cumulative exposure to sub-concussive head impacts was monitored for all participants using helmet-based accelerometers (gForce Tracker, GFT; Hardware version GFT3S ver4.0, Artaflex Inc., Markham, ON, Canada), which provided measures of linear acceleration (g) and rotational velocity (°/s), by impact location. The head impacts were categorized as ‘top’, ‘front’, ‘right’, ‘left’, and ‘back’. Sensors were mounted inside the helmet shell, to the left of the crown air bladder. ²⁴ Sensor data were collected during all games and practices (including training camp), using a trigger threshold of 15 g, in order to comply with existing literature validating the GFT system. ²⁴

Respiratory manipulations

Each session consisted of two functional MRI acquisitions collected separately with different breathing challenges: one hypercapnia (HC), and one hyperoxia (HO). The breathing manipulations consisted of a 2-min baseline period, followed by a block of gas inhalation (2-min for HC and 3-min for HO), and a 2-min recovery period (Figure 2(a)). A break of approximately 2 min was given to each subject between the protocols, during which additional imaging was acquired (M₀; see below). HC and HO blocks were induced using a feed-forward computer-controlled gas blending system (RA-MR™, Thornhill Research Inc., Toronto, ON), connected to a sequential gas delivery circuit. End-tidal pressures of carbon dioxide (P_ETCO₂) and end-tidal pressure of oxygen (O₂) (P_ETO₂) were sampled continuously via lines connecting the subject to the automated system (Figure 2(a)). Prior to entering the magnet bore, subjects were fitted with a facemask sealed using adhesive tape (Tegaderm, 3 M Health Care, St. Paul, MN, US). Subjects were then familiarized with the hypercapnic breathing paradigm using a condensed version of the protocol (1 min). OPEN IN VIEWER

During HC, P_ETCO₂ was targeted at 10 mmHg above the subject’s baseline using a hypercapnic gas mixture blended with medical (21:79% oxygen (O₂):N₂) and carbon dioxide (CO₂) enriched air (10:21:69% CO₂:O₂:N₂). P_ETO₂ was clamped at a baseline value of 110 mmHg. During HO, P_ETO₂ was targeted at 300 mmHg, while P_ETCO₂ was maintained at baseline values (Figure 2(a)).

MRI data acquisition

Imaging data were acquired using a Siemens 3.0 T Magnetom Tim Trio system with a 32-channel receiver head coil. A whole-brain T₁-weighted MP-RAGE (magnetization prepared rapid acquisition gradient echo) image was acquired for segmentation and registration using the following parameters: TR = 1760 ms, TE = 2.2 ms, time of inversion (TI) = 900 ms, voxel size = 1 mm isotropic, field of view (FOV) = 256 × 256 mm, flip angle = 9°, receiver bandwidth = 200 Hz/pixel, scan duration = 7 min and 32 s.

BOLD and CBF data were acquired simultaneously throughout the breathing manipulations using a dual-echo pseudo-continuous arterial spin labeling sequence (pCASL) ²⁵ with echo planar imaging (EPI) readout: TR = 4000 ms, TE₁/TE₂ = 10/30 ms, FOV = 250 ×250 mm, flip angle = 90°, voxel size = 3.9 mm isotropic, slice gap = 0.773 mm, label offset = 100 mm, receiver bandwidth = 2604 Hz/pixel, and EPI factor = 64. The label duration was set to 1665 ms, while the post-labeling delay (PLD) ranged from 1000 to 2291 ms, based on the multislice single-shot EPI imaging readout, which caused a different PLD in each slice acquired (average PLD = 1646 ms). This correction was accounted for in the computation of CBF maps (see below). Twenty-five axial slices were acquired on a 64 × 64 matrix (7/8 partial Fourier) in ascending order (i.e. bottom to top) for whole brain coverage using parallel imaging (GRAPPA acceleration factor = 2). For quantification of CBF, a tissue equilibrium magnetization map (M₀) was acquired using the same pCASL parameters, with a longer TR (15,000 ms) and no spin labelling. EPI images were calibrated using a pre-scan normalized image acquired using the body transmit/receiver coil. This step served to correct for inhomogeneities in the receive sensitivity of the 32-channel head coil.

Data preprocessing

Data from this study were preprocessed using a combination of FSL, ²⁶ AFNI, ²⁷ and in-house designed Matlab (MATLAB 2018b, The MathWorks, Inc., Natick, Massachusetts, USA) scripts. Prior to reconstruction of each signal, non-brain tissues were removed using the brain extraction tool (BET), followed by automated motion correction (FSL MCFLIRT).^28,29 FSL’s epi_reg tool ²⁶ was then used to co-register and align each subject’s HO and HC in native space, while combining topup^30,31 and boundary-based registration ³² to improve alignment between functional scans. Following automatic 3D segmentation of the structural T₁ image (FAST), ³³ the individual’s extracted GM mask was re-sampled linearly into low-resolution native space using a reversed rigid-body transformation matrix (6 dof; FLIRT) ²⁸ that aligned the EPI image to the structural anatomical scan.

BOLD images were isolated from the corrected pCASL series using a surround averaging of the second echo series (TE = 30 ms). ³⁴ The BOLD data were then spatially smoothed using a Gaussian kernel of 8 mm with SUSAN, ²⁶ and high-pass filtered ²⁶ to address signal drift during the EPI acquisition. CBF data were reconstructed using a linear surround subtraction of the first-echo ASL image (TE = 10 ms), between the adjacent tag and control frames and spatially smoothed (8 mm full width at half maximum Gaussian kernel). To account for confounding effects of motion artefacts on ASL signal, and possible over-estimation of CBF, the BOLD and perfusion timeseries from each breathing paradigm were despiked (3dDespike, AFNI) using an upper range of 2.5 standard deviations. The perfusion data were then converted into absolute physiological units (mL/100 g tissue/min) using FSL’s oxford_asl toolbox, ³⁵ with partial volume correction, ³⁶ in order to reduce the partial volume effects from neighboring voxels containing a mixture of GM, white-matter and cerebrospinal fluid, as well as the possible confounding effects of changes in tissue integrity over time. The ‘t2star’ option was also selected to correct for T₂* effects in the computation of CBF maps. The differences in PLD between each slice (i) due to the two-dimensional EPI readout were accounted using PLD_i = 1000 ms + (s_T) × (i − 1), where s_T represents the slice time in seconds (s_T = 53.8 ms). For images acquired during baseline and HO breathing, the inversion efficiency (α_inv) was set to 0.84. During HC, α_inv was set to 0.80 to account for the higher blood flow velocity in the feeding arteries, as a response to the higher arterial CO₂ content. ³⁷ The vasoconstrictive effects of HO on CBF ³⁸ were obtained from the dual-echo imaging data and incorporated into the calibration model. To account for the reduction in arterial blood T₁ (T_1b) during HO, due to higher plasma concentration of paramagnetic O₂,^39,40 the T_1b of blood was adjusted on a per-volume basis during CBF quantification based on the method proposed in Germuska et al. (Appendix A.1.2) ⁴¹ and the linear relationship between R₁ and PaO₂. ⁴²

Grey-matter volume maps were computed using steps extracted from the voxel-based morphometry pipeline in FSL. ³¹ Briefly, segmented grey-matter tissues from the T₁ anatomical scans were spatially transformed to the GM ICBM-152 template using affine (12 dof), ²⁸ and non-linear warp-fields (FNIRT). ⁴³ The spatially normalized maps were then used to create a symmetric study-specific grey-matter template, in order to avoid biases in the analysis based on morphological differences in tissues across subjects. Individual grey-matter maps were finally re-registered non-linearly to the grey-matter template and multiplied by the Jacobian warp field (‘modulation step’), followed by smoothing. These maps were then used for the volumetric statistical analysis of the grey-matter, over time (see below), based on exposure.

Computation of hemodynamic parameter maps

The first 5 volumes of each pCASL timeseries were removed to allow the respiratory pattern to reach steady-state. Volumes acquired during baseline and stimulus periods were averaged (Figure 2) to create mean baseline, hypercapnic and hyperoxic images. Next, the BOLD and CBF signal changes associated with the step HC and HO (Figure 2(b) and (c)) were used to calculate the M parameter (the theoretical maximal percent change from complete removal of deoxyhemoglobin), on a voxel by voxel basis, described in the original BOLD calibration model ⁴⁴

M = Δ BOLD BOLD 0 1 - ( CBF CBF 0 ) α · ( [ dHb ] [ dHb ] 0 ) β

(1)

In equation (1), and hereafter, BOLD₀, and CBF₀ represent baseline values, whereas variables without subscripts represent fractional changes estimated from the HC and HO respiratory manipulation. These were calculated using the final 80 s period of the stimulus block (120 s), given that delays in the vascular response may affect the rapid transition period from baseline.^45–47 The Grubb coefficient (α), representing the non-linear coupling between flow and venous volume, ⁴⁸ was set to 0.18. ⁴⁹ The parameter β, which represents the non-linear relationship between changes in deoxy-hemoglobin concentration ([dHb]) and R₂*, was set to 1.3 (appropriate for field strength of 3.0 T). ⁵⁰

Baseline normalized [dHb] is a function of the changes in arterial oxygen content (C_aO₂), OEF₀, and fluctuations in flow throughout the breathing protocols via equation (2) ⁵¹

[ dHb ] [ dHb ] 0 = ( Ca 02 | 0 OEF 0 x 1 - Ca 02 | 0 x ( 1 - OEF 0 ) ) ( CBF 0 CBF ) + ( 1 - Ca 02 x 1 - Ca 02 | 0 x ( 1 - OEF 0 ) )

(2)

In equation (2), the O₂ carrying capacity of hemoglobin (ϕ) was set to 1.34 mL O₂/g (Hb), while the total concentration of hemoglobin ([Hb]; oxygenated and deoxygenated) was assumed to be 15 g (Hb)/dL blood. Resting (C_aO_2|0) and stimulus (C_aO₂) oxygen content were derived using the Severinghaus equation ⁵² to get arterial hemoglobin saturation (S_aO₂) and equation (3), which represents the sum of the oxygen bound to hemoglobin (first term) and dissolved in plasma (second term)

C a O 2 = ( ϕ · [ Hb ] · SaO 2 ) + ( PaO 2 × ɛ )

(3)

In equation (3), P_aO₂ was replaced by the mean P_ETO₂ for each block, and ɛ (0.0031 mL O₂/dL blood/mm Hg), which represents the solubility of O₂ in plasma. ⁵¹

M (equation (1)) and OEF₀ (equation (2)) were solved simultaneously on a voxel per voxel basis using the Gauthier method ⁵¹ to find a unique solution at the intersection of both functions, while integrating data from each breathing manipulation (see Figure 1 in Gauthier and Hoge ⁵¹ and Supplementary Figure 1(a) to (d)). Due to the low SNR in ASL, the low-resolution GM mask was used to limit the application of this method to GM voxels only, which are likely to provide better signal. Because flow changes in HO should be minimal,^46,53,54 the GM-averaged T₁-corrected CBF for baseline and HO were used as inputs in equation (2) for the hyperoxic challenge. CMRO₂ was assumed to be constant ( CMRO 2 / CMRO 2 | 0 = 1 ) in this model (iso-metabolism). ⁴⁹

From here, individual voxelwise CMRO₂|₀ maps corrected for partial volume effects were computed using

CMRO 2 | 0 = CaO 2 | 0 · CBF 0 · OEF 0 · [ 39 . 34 umol ml ]

(4)

Data analysis

Functional maps (e.g. CBF₀, OEF₀, and CMRO_2|0) were warped into 2 mm Montréal Neurological Institute (MNI) space using concatenated linear and non-linear transformations that registered native images to standard space, via the high-resolution anatomical scan. Prior to alignment onto the standard template, the individuals’ GM masks were filtered in low-resolution space using a T₂* map, which was derived from fitting the native multi-echo log-transformed signal with a linear regression ( T 2 * = - Δ TE / - TE log ( Signal TE 2 / Signal TE 1 ) log ( Signal TE 2 / Signal TE 1 ) ; Supplementary Figure 2(a)). ⁵⁵ This was done in order to eliminate voxels with baseline T₂* values below 30 ms, which are likely caused by susceptibility artefacts (Supplementary Figure 2(b)), and thus may bias the analysis of hemodynamic parameters over time. The subject’s filtered GM mask (thresholded at 40%) was finally resampled into 2 mm MNI space, for further analysis of the physiological parameters.

Prior to the statistical analysis, two exposure groups were determined based on the daily exposure data recorded from the helmet accelerometers, between the PRE and PTC time points. Following computation of the average count of threshold impacts per session, the median of the sample was used to split the 44 participants into two equally sized groups: ‘LOW’ (below the median) and ‘HIGH’ (above the median) exposure. The LOW group served as the control for this analysis, in order to account for possible variations in CBF₀, OEF₀, and CMRO_2|0 over time that may be independent of exposure to sub-concussive impacts.

Volumetric and corrected CBF₀ maps were assessed for significant changes over time using a two by three factors mixed ANOVA, in order to account for each time point (within-subject; PRE, PTC, POST), both exposure groups (between-subject; LOW, HIGH), and missing data throughout the season due to injuries. This analysis was restricted to filtered GM voxels only. Significance from the voxelwise analysis was corrected for family-wise error at P < 0.05 and a parameter-specific cluster size determined using Monte Carlo simulations (10,000 iterations) in AFNI’s 3dFWHMx (with spatial AutoCorrelation) and 3dClustSim. ⁵⁶ Significant clusters from the interaction between time and exposure were then used to extract regional volumetric and CBF₀ measurements and run post-hoc tests to understand the direction of the group differences based on time. Significant clusters from the CBF₀ analysis were also used to extract regional OEF₀ and corrected CMRO_2|0 values from each subject, in order to test for co-localized changes in metabolism, based on the effects of both exposure and time. Prior to the analysis of the regional values extracted from the clusters, the normality of each parameter’s distribution was tested using the Shapiro-Wilk W-test. If not normally distributed, parameter values for each subject were log transformed before statistical comparison in IBM SPSS statistics (version 24.0, SPSS Inc., Chicago, IL, USA). Baseline P_ETCO₂ values for each time point were used as a time-varying covariate variable when comparing CBF₀ and OEF₀.

Data-driven grouping for exposure based on helmet biometrics

The helmet accelerometer data were used to split subjects into a LOW and HIGH exposure group based on a median average frequency of impacts per session equal to 9.21 (Figure 3(a)). The LOW group (mean_age = 20.2 ± 1.2 years, mean_weight = 91.3 ±13.2 kg, mean_height = 184.8 ± 5.4 cm) consisted of eight defensive backs (DB), four defensive linemen (DL), one fullback (FB), one kicker (K), three linebackers (LB), one quarterback (QB), two safeties (S), and two wide-receivers (WR).The HIGH group (mean_age = 19.7 ± 1.6 years, mean_weigth = 96.3 ± 11.8 kg, mean_height =185.8 ± 5.7 cm) included two DBs, five DLs, seven LBs, one offensive lineman (OL), one FB, one running back (RB), one S, two tight-ends (TE), and two WRs. The maximum number of impacts per session was also compared between the groups to show that HIGH exposure players sustained greater high impact sessions, relative to the LOW group (average maximum impacts in single session for LOW = 23.3 ± 10.4 impacts and HIGH = 42.6 ± 21.9 impacts; P < 0.0001). There were no significant differences (P > 0.05) in age, height or weight between the two exposure groups (i.e., LOW vs. HIGH). OPEN IN VIEWER

Further statistical analysis of the helmet kinematic data based on the LOW vs. HIGH exposure groups showed that the HIGH exposure players also sustained greater average peak rotational velocity per session (LOW = 413.7 ± 54 °/s and HIGH = 456.5 ± 72 °/s; P = 0.018), although no differences in mean linear acceleration per session were documented between the exposure groups (LOW = 31.1 ± 3.8 g and HIGH =32.4 ± 3.5 g; P = 0.079). Finally, statistical analyses of the count of impacts per location also revealed that the HIGH exposure group sustained more impacts on all sides of the helmet per session recorded, compared to the LOW exposure players (Figure 3(b)).

Respiratory data

No significant differences in baseline P_ETCO₂ and P_ETO₂ were documented across the three time points, for each group, except small fluctuations (±3 mmHg) in P_ETO₂ for the LOW exposure players (Table 1). Additionally, subjects from each group were exposed to similar hypercapnic and hyperoxic gas challenges during the breathing manipulations, indicated by the constant magnitude of change in P_ETCO₂ and P_ETO₂ across all time points (Table 1).

OPEN IN VIEWER

Table 1.

End-tidal pressures of CO₂ and O₂ during baseline, hypercapnia, and hyperoxia in each exposure group.

	LOW			P	HIGH			P
	PRE	PTC	POST		PRE	PTC	POST
Hypercapnia
Baseline PETCO2 (mmHg)	40 (4)	40 (4)	39 (4)	0.051	41 (4)	40 (3)	40 (3)	0.936
Baseline PETO2 (mmHg)	110 (5)	111 (7)	108 (3)	0.101	108 (6)	106 (5)	107 (3)	0.556
ΔPETCO2 (mmHg)	9 (1)	9 (1)	9 (1)	0.871	8 (1)	8 (2)	9 (2)	0.169
ΔPETO2 (mmHg)	3 (3)	2 (6)	2 (2)	0.450	5 (3)	4 (4)	3 (3)	0.257
Hyperoxia
Baseline PETCO2 (mmHg)	40 (5)	40 (4)	39 (5)	0.248	40 (4)	40 (3)	41 (4)	0.395
Baseline PETO2 (mmHg)	111 (6)	110 (3)	108 (3)	0.041^	108 (7)	105 (7)	108 (4)	0.065
ΔPETCO2 (mmHg)	−1 (1)	0 (1)	0 (1)	0.682	0 (1)	−1 (1)	0 (1)	0.717
ΔPETO2 (mmHg)	240 (40)	248 (39)	231 (43)	0.097	239 (46)	238 (39)	230 (42)	0.846

Volumetric differences based on exposure and time

Overall, there were no significant interaction between exposure and time on the SCAT-3 symptom scale (P = 0.473) and severity scores (P = 0.128).

Voxelwise volumetric analysis of the grey-matter volume maps revealed a significant interaction between exposure and time over widespread regions of the brain (Figure 4). Five significant clusters (corrected for a minimum size of 148.6 voxels) were localized across both the left and right hemispheres, with specific focal points spread around the frontal, parietal, temporal, and occipital lobes (Table 2). OPEN IN VIEWER

OPEN IN VIEWER

Table 2.

Significant clusters from the voxelwise 2 × 3 mixed ANOVAs on the volumetric and cerebral blood flow data showing a significant interaction between time and exposure.

Cluster	Size (voxels)	Center of mass (mm) a			Location
		X	Y	Z
Volumetric data (grey-matter volume)
1	546	−43	−82	12	Left middle occipital gyrus
2	372	6	−79	24	Right cuneus and right Brodmann area 18
3	276	2	−25	40	Right cingulate gyrus
4	187	−66	−42	2	Left middle temporal gyrus and left Brodmann area 22
5	178	36	60	−2	Right superior frontal gyrus and right Brodmann area 10
Cerebral blood flow (with partial volume correction)
1	1417	−32	5	49	Left middle frontal gyrus and left Brodmann area 6
2	1069	26	−19	62	Right precentral gyrus and right Brodmann area 6
3	816	−34	−71	42	Left superior parietal lobule and left Brodmann areas 7 and 19

a

Coordinates were converted from Talairach space to the Montreal Neurological Institute template and identified using the 3dclust and whereami functions in AFNI. Clusters were thresholded at P < 0.05 and a minimum size of 149 and 797.5 voxels for grey-matter volume and cerebral blood flow, respectively.

Post-hoc mixed ANOVA testing for the effect of time within each group indicated a difference in grey-matter volume for the HIGH exposure players over time (P = 0.026). Further pairwise comparisons within this group showed that grey-matter volume was reduced (P = 0.009) at POST (0.53 ± 0.04 mm³), in comparison to the PRE-season time point (0.55 ± 0.03 mm³). No significant differences were documented in the LOW group (P = 0.226).

Exposure-specific changes in hemodynamic parameters over time

T2₀* filtering of the grey-matter maps removed approximately 4 ± 1% of voxels, on average, across subjects. These were more often located within the temporal and inferior frontal lobes.

Voxelwise analysis of CBF₀ maps corrected for partial volume effects showed a significant interaction for time and exposure (Figure 5(a)) in three large clusters (Table 2) spread over the frontal lobe and the superior region of the right and left parietal lobes. Post-hoc analysis of the regional perfusion measurements for the main effect of time showed that significant differences were driven by alterations in CBF₀ within the HIGH group (Figure 5(b); P = 0.001). More specifically, CBF₀ was decreased at the POST time point, compared to both the PRE (P = 0.002) and PTC (P < 0.0001) measurements. CBF₀ at PTC was not statistically different compared to the PRE-season values (P = 0.203) when corrected for the confounding effects of time-varying P_ETCO₂ (P < 0.0001). OPEN IN VIEWER

In order to assess the physiological mechanisms underlying changes in perfusion, regional OEF₀ and CMRO_2|0 measurements were also extracted from the voxels highlighted in Figure 5(a). No significant interaction between exposure and time was observed for changes in OEF₀, and thus no further post-hoc tests were conducted. This was true for statistical testing of OEF₀ with (P = 0.521) and without (P = 0.596) corrections for the confounding effect of P_ETCO₂ over time (P < 0.0001). Despite no change in OEF₀, significant fluctuations in regional CMRO_2|0 were identified across time points, based on exposure (Figure 5(b); P = 0.029). More specifically, CMRO_2|0 measurements were found to change significantly over time for the HIGH group (P = 0.020), with higher CMRO_2|0 values observed at the PTC time point, compared to both the PRE (P = 0.035) and POST (P = 0.008) values. No significant differences were recorded between the PRE and POST time point (P = 0.345). Lastly, no significant effect of time was documented within the LOW exposure group (Figure 5(b); P = 0.243).

Main findings

The present study is the first to combine volumetric and quantitative calibrated MRI, along with helmet accelerometers, to explore the effects of repeated sub-concussive impacts on brain structure and cerebral physiology. The main findings of this study are three-fold: (1) In the HIGH group, significant changes in grey-matter volume were observed following the season, after one month of no contact activity, suggesting that brain structure may be altered longitudinally following exposure to head impacts. (2) Decreased CBF₀ after the season, corrected for partial volume differences, was specific to players with greater daily exposure to sub-concussive impacts. (3) Significant increases in CMRO_2|0 mid-season were observed in the HIGH group, which, following one month post-season, returned to baseline values. Altogether, these findings support the hypothesis that exposure to repeated sub-concussive impacts to the head may be associated with sub-acute changes in brain structure and resting physiology, emphasizing the need to monitor the frequency of such impacts in football athletes throughout the season.

Time-specific fluctuations in structural integrity and resting cerebrovascular physiology

Currently, the literature on the structural and neurophysiological effects of exposure to sub-concussive impacts is limited. ⁵⁷ However, some studies have demonstrated that subtle alterations in hippocampal volume, ⁵⁸ cognitive function, ⁵ and biochemical and functional markers^59–62 can all be found in collision sport athletes, despite the lack of a diagnosed concussion. In this study, decreases in grey-matter volume were identified post-season in the HIGH group. This is in line with previous findings in soccer players,^63,64 suggesting that exposure to head impacts may be associated with structural differences longitudinally. Although some changes in hippocampal volume pre- and post-season were reported in Slobounov et al., ¹¹ these were not statistically significant, making findings from this study the first to be put forward in a cohort of collegiate football athletes and contribute to the literature on the effects of sub-concussive impacts on brain structure. Differences in grey-matter volume may indicate changes in the cellular integrity of the tissues, ⁶⁵ although additional animal work is required to better understand the relationship between repeated microtrauma and changes in tissue morphology. ⁷

In addition to structural changes in the grey-matter, changes in perfusion were also specific to players within the HIGH group, suggesting that sub-concussive head impacts may be associated with transient changes in resting physiology. Here, lower resting perfusion was observed at the recovery time point (POST), relative to both PRE and PTC time points, suggesting that changes in regional brain perfusion persisted at one month following the end of the season, and no exposure to contact. No significant changes in CBF₀ with respect to baseline were documented at the PTC time point, which differs from findings in Slobounov et al., ¹¹ who reported global and regional increases in resting perfusion following exposure to sub-concussive impacts. Differences in CBF₀ reported in Slobounov et al., ¹¹ compared to this study, may be confounded by the differences in the length of exposure between the baseline and follow-up MRIs, which was much shorter in this study (∼1 month), as opposed to an entire season in Slobounov et al. ¹¹ Their design also did not account for time-varying changes in end-tidal CO₂, which was found to covary significantly with changes in CBF₀, in this study.

CMRO_2|0 was increased in the HIGH group after onset of participation in contact activities, which was in line with our original hypothesis. These results suggest that football players in the HIGH group had greater resting metabolic demand in-season, compared to both baseline and post-season time points. Group-specific differences in CMRO_2|0 over time contribute evidence for the hypothesis that the frequency of sub-concussive impacts is related to changes in resting cerebral physiology, proposed by Slobounov et al. ¹¹ Higher regional CMRO_2|0 at the PTC time point may result from regional increases in cellular O₂ uptake due to inefficiencies in oxygen utilization, ⁶⁶ or, greater consumption of energy (in the form of ATP) in order to meet elevated metabolic demands in response to shear stresses, ¹⁴ or trauma from mechanical stretch of the neurons. ⁶⁷ In both scenarios, the observed increase in CMRO_2|0 would ensure that the brain is supplied with adequate energy pre-cursors, which may be directed towards re-establishing ionic homeostasis and membrane potential disturbed following cellular damage.^68–70 While elevated CMRO_2|0 measurements in the HIGH group returned towards baseline values one month following the end of the season, CBF₀ was significantly reduced post-season within the region of interest (Figure 5(b)). In healthy tissues, decreases in CBF₀ are typically counteracted with higher OEF₀, in order to maintain homeostatic balances between delivery, and consumption of oxygen. Here, in the HIGH group, persistent OEF₀ values POST-season, in line with measurements from PRE and PTC, suggest that decreases in flow (or maintained OEF₀ values) may be a physiological response from homeostatic mechanisms attempting to balance cerebral metabolic demands.

Using helmet biometrics as an index for subject-specific differences in exposure to head impacts

Knowledge about the differences in exposure to head impacts across football positions^71,72 has been used to classify players and explore the relationship between microtrauma and acute or long-term changes in brain biomarkers.^6,73 Although this is common practice, the distribution of positions between the LOW and HIGH group in this study suggests that careful consideration should be made when assigning athletes to different groups based on their predicted exposure to head impacts. It is likely that certain positions such as OLs, DLs, and LBs may often be assigned to the HIGH exposure group. ⁷² However, there also seems to be a certain degree of unpredictability and overlap in positions like WRs, DBs, and even DLs, who vary widely in terms of their exposure to head impacts per session, likely due to individual differences in playing style and technique. Additional factors like differences in aggressivity, ⁷⁴ and/or starting status, ¹ can also play a critical role in determining how often players sustain impacts to the head, during a certain time period. Thus, moving forward, neuroimaging studies of head impacts should strongly consider using helmet telemetry in order to properly identify hitting profiles within each athlete, and adequately separate participants based on their own individualized behavior, and not the expectation for a position, or body size.

Limitations and future considerations

During the calibration, [Hb] was assumed to remain constant across all time points throughout the season. However, possible variations in [Hb] over time, or between the groups, based on changes in hematocrit levels due to fitness ⁷⁵ or mental stress^76,77 for instance, may also have contributed to the observed changes in resting cerebral metabolism. This effect was likely minimized in this design given that group-differences documented over time were extracted based on exposure to head impacts in players from varying positions, included in each experimental group (LOW vs. HIGH; see section ‘Participants’ section). Little to no differences in CMRO_2|0 (or CBF₀) were observed in players who sustained lower exposure to sub-concussive impacts during the timeline studied. Furthermore, because in-season training is heavily focused on strength and fitness maintenance built during the off-season, possible fluctuations in [Hb] over time may have been limited. Future designs should, however, aim to collect a sub-sample of blood in athletes studied, over time, in order to provide more insight about possible variations in [Hb] during the season and improve the calibration approach.

Other modeling parameters such as α and β were chosen based on existing literature, ⁵¹ although more recent methods have presented novel ways to estimate such parameters. ⁷⁸ These methods may improve the estimation of OEF₀, although they also present with methodological limitations and still require additional validation across different clinical populations. In this design, the α constant was used to estimate changes in cerebral blood volume (CBV) indirectly via changes in CBF during the gas manipulations, based on the power–law relationship. ⁴⁸ The lack of direct CBV measurements limits our ability to estimate regional CMRO₂ using the biophysical model, as α may vary depending on the experimental conditions. ⁷⁹ Future designs may consider the use of vascular space occupancy (VASO) ⁸⁰ in order to determine CBV experimentally, ⁸¹ and possibly improve calibrated CMRO₂ measurements. ⁸²

In this study, a nominal PLD of 1000 ms was used to acquire the ASL signal, which was adjusted on a slice-by-slice basis due to the 2D EPI readout. ASL measurements of CBF are dependent on the selection of the PDL, which inherently relates to the longitudinal relaxation time of blood, ⁸³ and the arterial transit time (ATT) from the labeling slice to the brain tissues.^83,84 ATT may vary across different regions of the brain, or between participants, which should be acknowledged when using ASL perfusion-based imaging in clinical population. No background suppression was employed in this imaging acquisition, although future design should also consider implementing such approach in order to improve SNR of the ASL signal. ⁸⁵

The effects of HC and HO on CMRO₂ were assumed to be negligible in this study, although this remains a topic of debate.^86–89 Changes in CMRO₂ during HC would violate an important assumption of the calibration model, which is that HC is assumed to be a purely vascular stimulus. Changes in CMRO₂ would propagate into erroneous estimation of the M parameter, and thus, OEF₀. However, findings from the original calibration method suggests that potential changes in CMRO₂ during HC are likely to be minimal. ⁵¹ This is further supported by the agreement between calibrated fMRI results from the studies reported above and the normative values published in the PET literature.^90–92

In comparison to work using PET, ⁹³ estimates of OEF₀ and CMRO_2|0 using calibrated fMRI have higher intra- and inter-subject variability,^55,94 likely due to the error propagation of the calibration. Reasons for this variability include the magnitude of the respiratory change in end-tidal measures, ⁹⁵ along with variations across subjects, and between time points, in the vasoactive and vasoconstrictive effects of HC and HO on CBF, respectively. In earlier methods,^44,51,54,96 a T₁-corrected percent change in CBF was assumed based on the entire group’s average decrease in flow, during HO. In this study, we estimated changes in arterial blood T₁ during HO, and corrected CBF on an individual basis, which likely helped minimize potential intra-subject variation across time points. Fluctuations in P_ETCO₂ between subjects and time points were also accounted for in our statistical model. Lastly, the effects of possible differences in OEF₀ and CMRO_2|0 within individuals, independent of exposure to head impacts, were also accounted for using the LOW group, which showed no significant variations over time, within the region of interest.