Scale‐free functional brain dynamics during recovery from sport‐related concussion

Abstract Studies using blood‐oxygenation‐level‐dependent functional magnetic resonance imaging (BOLD fMRI) have characterized how the resting brain is affected by concussion. The literature to date, however, has largely focused on measuring changes in the spatial organization of functional brain networks. In the present study, changes in the temporal dynamics of BOLD signals are examined throughout concussion recovery using scaling (or fractal) analysis. Imaging data were collected for 228 university‐level athletes, 61 with concussion and 167 athletic controls. Concussed athletes were scanned at the acute phase of injury (1–7 days postinjury), the subacute phase (8–14 days postinjury), medical clearance to return to sport (RTS), 1 month post‐RTS and 1 year post‐RTS. The wavelet leader multifractal approach was used to assess scaling (c1) and multifractal (c2) behavior. Significant longitudinal changes were identified for c1, which was lowest at acute injury, became significantly elevated at RTS, and returned near control levels by 1 year post‐RTS. No longitudinal changes were identified for c2. Secondary analyses showed that clinical measures of acute symptom severity and time to RTP were related to longitudinal changes in c1. Athletes with both higher symptoms and prolonged recovery had elevated c1 values at RTS, while athletes with higher symptoms but rapid recovery had reduced c1 at acute injury. This study provides the first evidence for long‐term recovery of BOLD scale‐free brain dynamics after a concussion.


2011
). In the context of sport, the diagnosis and management of concussion are mainly based on symptom assessments, along with brief cognitive and balance testing, with the clinical determination of return to sport (RTS) based on symptom resolution following a graded exercise protocol . Although the natural history of clinical recovery is well described, brain recovery following concussion is less understood (McCrea et al., 2017), particularly with respect to time of RTS.
As a mild form of traumatic brain injury (TBI), concussion is rarely associated with gross neuroradiological findings. Instead, injuries at the microscopic level, undetectable in individual patients using standard diagnostic imaging, have a measurable impact on brain function, including disturbances in neurometabolic activity and the regulation of cerebral blood flow (Giza & Hovda, 2014). Blood-oxygenation-level dependent functional magnetic resonance imaging (BOLD fMRI) has been used to detect changes in brain function among concussed individuals, during cognitive tasks and at rest . Studies of the resting brain have typically focused on functional connectivity, which measures the synchrony of spontaneous BOLD signal fluctuations between different brain regions (Sporns, 2014). Functional connectivity has been studied in concussed cohorts, during the symptomatic phase of injury and after RTS (Churchill et al., 2017b;Johnson et al., 2012;Zhang et al., 2012;Zhu et al., 2015), showing significant concussion-related disturbances that persist beyond medical clearance.
Connectivity-based methods provide an incomplete picture of the functional brain changes that occur after a concussion, however.
These methods do not directly characterize the temporal dynamics of spontaneous BOLD signal fluctuations or provide information about how these dynamics are affected by injury. The temporal dynamics of BOLD signals have been well characterized in the uninjured resting brain and show evidence of scale free (or fractal) behavior, such that no specific timescale plays a dominant role (He, 2014). In the same way that geometric fractals are "self-similar," with each magnified part resembling a smaller-sized copy of the whole (Mandelbrot, 1983), a scale-free time series signal x(t) is statistically indistinguishable from a dilated and rescaled version of itself, that is, x(t)~a −H x(at), for all values of a > 0. In this expression, the Hurst exponent H parameterizes scaling behavior, with higher values indicating a more scale-free signal. In practice, the scaling behavior of functional brain signals is often defined in terms of a power-law relationship between frequency f and the power spectral density (PSD) P x (f ) / |f| −β , where β = 2H − 1 (Eke et al., 2000). More recently, sophisticated multifractal formalisms have extended beyond monofractal models in which a single scaling behavior H is modeled, instead describing signals in terms of a spectrum of time-varying scaling exponents h, thereby providing a richer description of BOLD dynamics (Ciuciu, Varoquaux, Abry, Sadaghiani, & Kleinschmidt, 2012;Shimizu, Barth, Windischberger, Moser, & Thurner, 2004; Wink, Bullmore, . Although it has long been known that the brain exhibits scale-free dynamics at multiple different levels of organization (Werner, 2010), the initial BOLD fMRI studies of scaling behavior considered it to be a product of measurement noise and/or physiological fluctuations unrelated to brain function . Subsequent studies challenged this assumption though, showing that scale-free processes in electrophysiology are linked to BOLD signal fluctuations (Van de Ville, Britz, & Michel, 2010). Moreover, BOLD scaling encodes important information, as scaling behaviors vary with cognitive state and systematic differences have been observed in brain areas serving different cognitive roles (Barnes, Bullmore, & Suckling, 2009;Ciuciu et al., 2012;He, 2011). Potentially more relevant to concussion, the suppression of BOLD scaling in the brain is associated with greater cognitive effort (Churchill et al., 2016), trait anxiety (Tolkunov, Rubin, & Mujica-Parodi, 2010), and distress during adverse life events (Churchill et al., 2015;Tolkunov et al., 2010). Reduced BOLD scaling is also associated with slower reaction time during tasks . In general, these findings point toward the suppression of BOLD scaling as an indicator of a more taxed, less adaptive brain. This is consistent with the literature on scaling in biological systems (Goldberger et al., 2002), which finds that reduced scaling, which reflects a loss of system complexity, is a hallmark of functional impairment.
Given the literature evidence that suppression of BOLD scaling is a marker of brain dysfunction, this study examined whether concussion, which involves impairments in brain function and behavior including cognitive difficulty and emotional dysregulation, had similar effects on BOLD scaling. It was hypothesized that for concussed athletes: (a) BOLD scaling will show continual increase from acute injury up to 1 month post-RTS, reflecting brain recovery that lasts beyond medical clearance; (b) for brain areas showing longitudinal change, BOLD scaling will be suppressed relative to uninjured controls at acute injury; and (c) longitudinal recovery will depend on clinical covariates of acute symptom severity and time to RTS. These hypotheses were evaluated using resting-state BOLD fMRI data acquired from a sample of university athletes with sport-related concussion who were followed longitudinally from the acute phase of injury to 1 year post-RTS. In addition, a large normative cohort of athletic controls was also evaluated. Scaling analyses were performed using the wavelet leader multifractal (WLM) formalism, which offers superior performance over standard monofractal techniques (Lashermes, Jaffard, & Abry, 2005;Wendt, Abry, & Jaffard, 2007).

| Study participants
Sixty-one concussed athletes were recruited from university-level sport teams at a single institution (including volleyball, hockey, soccer, football, rugby, basketball, lacrosse, and water polo; see Supplementary Table S1 for athlete numbers by sport) through the academic sport medicine clinic, following a concussion diagnosis. Diagnosis was determined by a staff physician following a sustained direct or indirect contact to the head with signs and/or symptoms as per the Concussion in Sport Group guidelines . Imaging evalua- (45/61), and 1YR (32/61). Attrition was not significantly related to demographic variables (age, sex, concussion history) or clinical variables (symptom severity, time to RTS) examined in this study, based on Spearman correlations, at a false discovery rate (FDR) of 0.05 across time points.
As a control group, one hundred and sixty-seven athletes were also consecutively recruited and imaged at the start of their competitive season. All athletes in the study completed baseline assessments with the Sport Concussion Assessment Tool (SCAT) (Echemendia et al., 2017;Guskiewicz et al., 2013)
Resting-state fMRI data were acquired via multislice T2*weighted echo planar imaging (TE/TR = 30/2,000 ms, FA = 70 , 32 oblique-axial slices with FOV = 200 × 200 mm 2 , 64 × 64 matrix, 4.0 mm slice thickness with 0.5 mm gap, 3.125 × 3.125 mm 2 in-plane resolution, BW = 2,298 Hz/px), producing a time series of 195 images at each slice location. During fMRI, athletes were instructed to lie still with their eyes closed and not to focus on anything in particular.
Processing and analysis were performed using the Analysis of Functional Neuroimages (AFNI) package (afni.nimh.nih.gov), fMRIB Software Library (FSL; https://fsl.fmrib.ox.ac.uk), and customized algorithms developed in the laboratory. After discarding the first four volumes to allow the fMRI signal to reach equilibrium, the processing included rigid-body motion correction (AFNI 3dvolreg), removal of outlier scan volumes using the SPIKECOR algorithm (nitrc.org/projects/ spikecor), slice-timing correction (AFNI 3dTshift), spatial smoothing with a 6 mm full width at half maximum (FWHM) isotropic 3D Gaussian kernel (AFNI 3dmerge) and regression of motion parameters and linear-quadratic trends as nuisance covariates. For motion parameter regression, principal component analysis was performed on the six rigid-body movement parameters (consistently accounting for >85% of variance), and the first two PCs were used as nuisance regressors.
To control for physiological noise, the data-driven PHYCAA+ algorithm (nitrc.org/projects/phycaa_plus) was used to downweight areas with nonneural signal, followed by further regression of signal originating from WM and cerebrospinal fluid (CSF). The WM and CSF regressions were performed after spatial normalization, described in the paragraph below.
To perform group-level analyses, the fMRI data were coregistered to a common anatomical template using the FMRIB Software Library (FSL) package (https://fsl.fmrib.ox.ac.uk). The FSL flirt algorithm was used to compute the rigid-body transform of the mean functional volume for each athlete to their T1-weighted anatomical image, along with the 12-parameter affine transformation of the T1 image for each athlete to the MNI152 template. The net transform was applied to the functional imaging data, which was resampled at 3 × 3 × 3 mm 3 resolution. To remove WM and CSF signal, subject T1-weighted images were segmented and coregistered to the MNI152 template using the fslvbm protocol (fsl.fmrib.ox.ac.uk/ fsl/fslwiki/FSLVBM), which used fast to obtain partial volume segmentation maps of gray matter (GM), WM, and CSF, followed by iterative applications of affine registration algorithm flirt and nonlinear registration algorithm fnirt, to obtain a symmetric, studyspecific mean GM tissue template. The spatial transforms were subsequently used to obtain mean WM and CSF tissue templates, resampled into 3 × 3 × 3 mm 3 resolution and a 6 mm FWHM isotropic 3D Gaussian smoothing kernel was applied. For WM, the brain mask p(WM) ≥ P 95% (WM) was obtained (i.e., voxels within the distribution 95th percentile) and a single spatial erosion performed (3 × 3 kernel, in-plane). Two mean seed time series were obtained by separately averaging over cerebral WM voxels and averaging over brainstem WM (as their time courses were substantially different). For CSF, the brain mask p(CSF) ≥ P 95% (CSF) was obtained and manually edited into two separate masks of the lateral ventricles. Two mean seed time series were obtained by separately averaging over these two ventricular regions. The four physiological time series were then regressed from each voxel, for all study participants.
To reduce computational burden and improve the stability of regional BOLD measures, the data were then parcellated using the Brainnetome Atlas (BNA), which a connectivity-based atlas that subdivides the brain into 210 cortical regions and 36 subcortical regions (Fan et al., 2016). Within each parcel, a mean seed time series of interest was obtained, for a cubic region of interest (ROI) of 3 × 3 × 3 voxels, placed at the parcel center of mass. Seed ROIs of uniform size were created for all parcels, to avoid potential unequal smoothing effects across caused by averaging over parcels of different size. This procedure generated a set of 246 BOLD time series per participant for subsequent analyses.

| Clinical and demographic data
Participant demographics are reported in Table 1, including age, sex, and prior concussion history, along with time to RTS for concussed athletes, defined as the number of days from concussion event to symptom resolution following a graded exertional protocol (McCrory et al., 2017). From the SCAT, a symptom severity score was obtained by summing across the 22-item symptom scale, with each item receiving a 7-point Likert scale rating. A total symptoms score was also obtained by counting the total number of symptoms with nonzero ratings. In addition, brief cognitive testing scores were reported, including orientation, immediate memory, concentration, and delayed memory, along with total scores for the modified balance error scoring system. All scores were tested at acute injury and RTS for a significant difference relative to baseline, via nonparametric Wilcoxon paired-measures tests (two-tailed). As per evolving clinical guidelines, the first 44/61 concussed athletes and 68/167 controls in this study were evaluated with SCAT3, while the remainder were evaluated using SCAT5. The immediate memory and delayed memory tests were changed in SCAT5 (from 15 to 30 items and from 5 to 10 items, respectively). For these subtests, statistics are only reported for SCAT3 data, as this represents the larger sample of concussed athletes. For all other subtests, statistics are based on the complete dataset, after verifying that there were no significant within-cohort differences in SCAT3 and SCAT5 scores based on two-sample Wilcoxon tests (p ≥ .288, for all measures).

| Measuring BOLD scaling
The goal of the analyses in this section was to summarize the dynamics of the spontaneous, arrhythmic fluctuations of restingstate BOLD time series x(t). This was attained for each of the 246 parcel ROIs identified for a given participant (concussed and control) and imaging session. One approach involves analyzing signal variance at different frequencies f using a PSD estimator. However, both the BOLD signal and underlying neural activity have a broad spectral distribution, where power declines smoothly with increasing frequency (Fox, Snyder, Vincent, & Raichle, 2007); hence, no specific set of frequencies or timescales can be singled out for analysis. For this reason, scaling analysis is an appealing alternative; instead of analyzing specific frequencies, this approach characterizes power-law scaling behaviors that relate the dynamics of BOLD signals at different timescales. By convention, x(t) is deemed scale invariant if the following relation holds for a wide range of frequencies: This has nontrivial implications, as the relative spectral power at a given frequency f is then entirely dictated by a single scaling exponent β. This parameter is also related to the Hurst exponent by β = 2H − 1 for a fractional Brownian motion model (Eke et al., 2000), where H is considered an index of temporal dependence. A value of 0 < H < 0.5 indicates short-range dependency, where a high value for x(t) tends to be followed by a low value for x(t + 1) and vice-versa; H = 0.5 denotes T A B L E 1 Demographic and clinical data for athletes with concussion and controls. Clinical scores of total symptoms and symptom severity are summarized by the median [Q1, Q3]. For tests of Immediate Memory and Delayed Memory, denoted by a "*," statistics are based on a reduced sample of 44/61 concussed athletes and 68/167 controls, due to changes in scoring guidelines between SCAT3 and SCAT5. Significant differences in scores at acute injury, relative to baseline, are noted with "**"

Control Concussion
Age ( an uncorrelated process; and 0.5 < H < 1.0 indicates long-range dependency, where a high (or low) value for x(t) tends to be followed by a high (or low) value for x(t + 1). In practise, PSD-based scaling analysis typically involves estimating β from the slope of log(P x (f )) versus log(f ) scatterplots, calculated using standard linear regression techniques.
Although intuitive, PSD-based scaling analyses suffer from a significant limitation, namely the inability to distinguish true scaling behavior from deterministic trends in the data. To address this issue, scaling may be evaluated using wavelets (Abry, Gonçalvés, & Flandrin, 1995;Abry & Veitch, 1998;Veitch & Abry, 1999). The wavelet transform employs translated and dilated versions of a basis function Ψ 0 (t), which is compact in both time and frequency domains, to analyze x(t) at different delays and timescales. This trends in x(t) (Daubechies, 1992). In the continuous time domain, the wavelet coefficient d x (a, k) measures signal at timescale a and delay k by computing the inner product d In discretely sampled real-world data, the analogous computation uses the discrete wavelet transform (DWT), which analyzes dyadic scales a = 2 j and delay intervals k = n2 j , for nonnegative integer j and n. The mean squared wavelet energy at a given timescale is then calculated The wavelet energy has a relationship with dyadic scale that is described by the following equation: The DWT approach shows reduced modeling error compared to PSD estimators (Abry et al., 1995;Abry & Veitch, 1998) and is thus the preferred technique for scaling analyses. As with PSD-based methods, scaling analysis typically involves estimating β from the slope of log(|d x (2 j )| 2 ) versus j scatterplots.
One additional concern is that the DWT approach assumes the data of interest to be Gaussian, with scaling of x(t) fully defined by its variance (i.e., its second-order moment) across timescales. We may extend the approach to encompass scaling behaviors for a range of moments q using the WLM formalism (Lashermes et al., 2005;. To ensure robust performance, wavelet coefficients d x (a, k) are replaced with wavelet leaders L x (a, k), which are calculated as the largest coefficient value | d x (a 0 , k 0 )| in a neighborhood of k, for all a 0 ≤ a. The WLM approach then describes a more general expression relating dyadic timescale j to wavelet power: Instead of a fixed β value, wavelet power scaling is now expressed in terms of a characteristic function τ(q) that varies with statistical moment q. This function is typically expressed as a Legendre polynomial expansion τ(q) ≈ P p c p (q p /p!), where the set of logcumulants c p summarize scaling behavior of x(t).
The function τ(q) is linked to the concept of multifractality, in which x(t) exhibits time-varying fractal features characterized by their Hölder exponent values h(t), which describe power-law scaling behavior of x(t) in the neighborhood of each time point t . This collection of scaling exponents is characterized by a "singularity spectrum" D(h), which measures the Hausdorff dimension at each value h, which is the "size" of the set of time points t i such that h(t i ) = h. In practise, D(h) is approximated using the Legendre transform of τ(q) and therefore the logcumulants c p reflect multifractal scaling behavior.
The present work examined the parameterization τ(q) = c 1 q + c 2 q 2 /2, as higher order coefficients were found to be nonsignificant under bootstrap testing . The log-cumulant c 1 is equivalent to the peak amplitude of D(h), that is, the most frequentlyoccurring scaling parameter h, which is similar to monofractal H . The log-cumulant c 2 determines the width of D(h), and can therefore be thought of as indexing the degree of multifractality (Ciuciu et al., 2012). Multifractal analysis was performed using the WLBMF toolbox (https://www.irit.fr/~Herwig.Wendt/software.html; ). The WLM estimation was performed using Daubechies' wavelets with N = 3 vanishing moments and q values ranging −10 to 10. Robust results were identified for dyadic wavelets over a three octaves range of j = [2, 4], which corresponds to a frequency range of [0.016, 0.125] Hz.

| Analysis of healthy controls
To characterize the spatial distribution of scaling behavior, the c 1 and c 2 log-cumulant values were obtained for all athletic controls, at each of the 246 ROIs. Bootstrap 1-sample tests were then conducted at each ROI to evaluate scaling behavior: c 1 was tested against the null c 1 ≤ 0.5 (1-tailed) and c 2 was tested against the null c 2 = 0 (2-tailed). Significant ROIs were then identified after adjusting for multiple comparisons at an FDR of 0.05. For ROIs that were determined to be significant, maps of the group mean and standard error values were displayed.  Given the demographic effects seen in the control cohort (see Section 2.5), a GLM was used to evaluate the effects of concussion with covariates adjusting for age, sex and history of concussion. This was done in a bootstrap resampling framework (1,000 iterations) to obtain coefficient b with a 95%CI, BSR, and p-value for each imaging session, with significant imaging sessions identified at an FDR threshold of 0.05. To mitigate potential bias and loss of efficiency due to missing data, bootstrap GLM analysis was combined with multiple imputation using the "Boot MI" approach of (Schomaker & Heumann, 2018): bootstrap samples were drawn from the full dataset (including missing data) and for each sample, imputation was done M = 10 times to generate 10 coefficient estimates, which were averaged to obtain a point estimate. The set of coefficient point estimates were treated as a conventional bootstrap empirical distribution, from which summary statistics were calculated. Imputation was done using the fitted LMM to generate simulated log-cumulant values. To evaluate the impact of imputation, results of MI analyses were also compared to unimputed bootstrap parameter estimates in Table S1.

| Effects of clinical covariates
A secondary set of analyses, performed within the concussed cohort, tested for effects of clinical covariates on log-cumulant values over the course of recovery. This included (a) total symptom severity at acute injury and (b) days to RTS. Given the high correlation between these variables (see Section 3.1), a pair of orthogonal composite scores (CS) were defined. After the two variables were renormalized via the inverse empirical distribution function and mean centered, composite score one (CS1) was computed as the average of the two variables (symptoms + days to RTP), which quantified overall clinical outcome. Composite score two (CS2) was computed as the difference (symptoms − days to RTP), which quantified discrepancy between the two measures of concussion outcome (i.e., a positive score denoting high symptom burden but rapid recovery, and a negative score denoting the converse). indicating that individuals with greater initial symptoms also tended to have a more prolonged recovery.  versus log(f ) and a steeper positive slope is evident in the plot of log 2 (d x (2 j )) versus j used for DWT estimation. Similarly, the plots of multifractal spectrum D(h) versus h obtained using WLM techniques show that the spectrum peak (i.e., c 1 ) of the precuneus is shifted to the right, denoting increased scaling behavior relative to the putamen. In addition, these spectra show relatively broad curves, F I G U R E 1 Illustration of BOLD scaling behavior, comparing example ROIs in the right precuneus (blue) and right putamen (red) for healthy control athletes. The GLM analyses of controls revealed significant effects of demographics on c 1 , with ROIs summarized in Table 2, whereas no significant effects on c 2 were identified. Among control athletes, greater age was significantly associated with increased c 1 in somatosensory and motor areas. Averaging over significant brain regions, a b coefficient value of 0.031 was obtained (95%CI: 0.020, 0.041; BSR = 6.18; p < .001). Sex was associated with more spatially extensive differences, as female athletes had lower c 1 values compared to male athletes throughout the brain, including prefrontal, occipital, and subcortical regions. Averaging over significant regions, a b value of −0.111 was obtained (95%CI: −0.142, −0.077; BSR = −6.73; p < .001). No significant effects of history of concussion on c 1 were identified.

| Main effects of concussion
The LMM analyses of concussed athletes showed significant effects of imaging session on c 1 , with ROIs summarized in Table 3, whereas no significant effects on c 2 were identified. Significant changes in c 1 relative to ACU were predominantly seen at RTS, with a single ROI (left medioventral occipital cortex), remaining significant at 1MO and two additional ROIs (left superior frontal and inferior frontal gyri) emerging as significant at this time; no significant effects were F I G U R E 2 Depiction of scaling coefficient values c 1 and c 2 for healthy control athletes, including group mean and standard error of the mean (SE). Regions of interest (ROIs) are shown as square patches and the MNI z-axis coordinate is listed below each slice T A B L E 2 Summary of ROIs identified in control athletes that show significant effects on c 1 for demographic factors of age, sex or history of concussion (conc.hx.). Standardized effect sizes are reported in terms of BSRs. Significant effects are noted with "*" (FDR = 0.05 threshold). Brain regions are defined based on the BNA T A B L E 3 Summary of ROIs identified in concussed athletes that show significant longitudinal change in c 1 relative to ACU, as displayed in Figure 4. Standardized effect sizes are reported in terms of BSRs. Significant effects are noted with "*" (FDR = 0.05 threshold). Brain regions are defined based on the BNA identified at 1YR. Figure 3a depicts the brain areas of significant longitudinal change. The ROIs were predominantly within occipital brain areas, with a smaller set in temporal and thalamic areas. All significant ROIs had uniformly positive BSRs, denoting a longitudinal increase in c 1 from ACU to RTS. Figure Table 4 also reports the two-sample analysis results comparing concussed athletes to controls. The c 1 values for concussed athlete were lower than controls at ACU, but the effects were nonsignificant. However, the c 1 values for concussed athletes were significantly elevated relative to controls at RTS, before becoming nonsignificant at later imaging sessions.
T A B L E 4 longitudinal effects of concussion on c 1 , averaged over all significant ROIs in Table 3. (left) longitudinal analysis using an LMM to compare postacute imaging sessions (SUB, RTS, 1MO, 1YR) to ACU. (Right) Cross-sectional analysis using a GLM to compare concussed imaging sessions (ACU, SUB, RTS, 1MO, 1YR) to uninjured controls. Statistics include coefficients of effect b, 95% confidence intervals (95% CIs), BSRs, and p-values. Significant longitudinal change relative to ACU is noted with "*" and significant cross-sectional difference relative to controls is noted with "**" (FDR = 0.05 threshold)  greater CS2 scores had lower c 1 values in these brain regions at early injury. Figure 4d plots the distribution of participant c 1 values as a function of imaging session, averaged over significant ROIs. For athletes with CS2 < 0 (lower symptoms, more days to RTS) average c 1 was above the average control value at ACU but declined toward control values in later imaging sessions, whereas for athletes with CS2 > 0 (higher symptoms, fewer days to RTS) average c 1 was below control values at ACU but increased toward controls in later imaging sessions.

| DISCUSSION
There is growing evidence that scaling analysis of functional brain dynamics captures important information about cognition and brain health. However, to date, most of the BOLD fMRI research in this domain has focused on healthy adults and has examined how scaling is modulated by cognitive tasks and other sources of effort (Barnes et al., 2009;Churchill et al., 2016;Ciuciu et al., 2012;He, 2011). The present study is the first to examine the relationship between scaling behavior in resting-state fMRI and sport-related concussion, with reference to a large control group. This study is also the first to show that functional brain dynamics are significantly altered in the course of concussion recovery, with effects that are modified by clinical variables of acute symptom severity and time to RTS.
The initial analyses of uninjured control athletes identified strong BOLD scaling behavior (c 1 ) with significant multifractality (c 2 ) throughout the healthy resting brain.  (Ciuciu et al., 2012) also reported the highest scaling in networks encompassing dorsoparietal, visual, and parieto-cingulate regions. Given the large control cohort in this study and the supporting literature, we may therefore conclude that BOLD scaling shows consistent spatial organization in the brain, with the highest values for cortical regions implicated in higher level processing, for example, executive function, attention, and interoception (Behrmann, Geng, & Shomstein, 2004;Critchley, Wiens, Rotshtein, Öhman, & Dolan, 2004;Fuster, 2000) and the lowest values in subcortical regions.
Further analyses of the controls identified significant demographic effects on scaling (c 1 ) but not multifractality (c 2 ). Greater age was associated with increased c 1 , albeit for a spatially limited set of ROIs, likely due to the narrow age band being studied. Nevertheless, results are consistent with a large-scale study of resting electrophysiology and brain maturation (Smit et al., 2011), which reported increasing long-range dependence with age until a plateau at 25-30 years, with similar trends in functional brain topology (Dosenbach et al., 2010). The observed effects were mainly localized in somatosensory and motor areas. One explanation is that athletes are improving skills involving these domains during their training, leading to increased neural efficiency (Guo, Li, & Yu, 2017;Wei & Luo, 2010) and therefore more scale-free signal among older athletes. In addition, female athletes exhibited lower BOLD scaling than male athletes, within a distributed set of prefrontal, occipital, and subcortical regions. Previous studies have reported sex differences in functional brain topology (Tian, Wang, Yan, & He, 2011;Zuo et al., 2010); however, to our knowledge, sex differences in resting-state BOLD scaling dynamics have not been previously examined. History of concussion had no significant effects on BOLD scaling, which aligns with nonsignificant findings in studies of resting-state functional connectivity and taskbased activation for athletes with a history of concussion (Churchill, Hutchison, Leung, Graham, & Schweizer, 2017;Terry et al., 2012).
Similarly, among concussed athletes, c 1 in posterior brain regions did not show significant effects of age but did show effects of sex that were comparable to uninjured controls. In contrast with the controls, however, concussed athletes with a history of concussion had significantly reduced posterior c 1 , suggesting a cumulative effect of repeated brain injury on BOLD dynamics during concussion recovery. This is consistent with prior literature, in which a history of concussion was associated with greater disturbances in functional connectivity and neurometabolites during recovery from subsequent sportrelated concussion Vagnozzi et al., 2008).
The analyses of concussed athletes identified significant longitudinal changes in BOLD scaling (c 1 ), whereas multifractality (c 2 ) showed no evidence of longitudinal change. These findings partly support our first study hypothesis, as BOLD scaling values were lowest at acute injury and increased in later imaging sessions. The findings are also congruent with prior literature, in which reduced BOLD scaling is associated with greater cognitive burden (Barnes et al., 2009;Churchill et al., 2016;Ciuciu et al., 2012;He, 2011) and mental distress (Churchill et al., 2015;Tolkunov et al., 2010). Unlike conventional BOLD measures of task-related activity and functional connectivity, which may show either increases or decreases, scale-free dynamics consistently show suppression associated with these cognitive and psychological stressors. This has been observed as a common feature of many biological systems, where a reduction in scaling behavior signals a loss of system complexity (Goldberger et al., 2002) and thus a more constrained, less adaptive system. The present study findings suggest that acute concussion, which is associated with impairments across multiple domains, including cognition (Echemendia, Putukian, Mackin, Julian, & Shoss, 2001) and emotion regulation (McCrory et al., 2017), may similarly involve disruptions in the scalefree dynamics of underlying brain function. However, our second study hypothesis was not supported, for although the c 1 values of concussed athletes were low relative to controls, the difference was nonsignificant. This is consistent with previous literature which has reported heterogeneous effects of concussion on resting brain function in the first week post-injury (Churchill et al., 2017a). These findings also reinforce that the effects of concussion on brain function are modest, compared to more severe forms of TBI (McAllister, Flashman, McDonald, & Saykin, 2006).
Intriguingly, post-RTS findings for concussed athletes differ from those predicted by our first study hypothesis. Instead of a continual increase in c 1 over time, the values for concussed athletes peaked at RTS, followed by a decline up to 1 year post-RTS. The c 1 values of concussed athletes at RTS were also significantly elevated compared to uninjured controls. These findings suggest that persistent concussion-related disturbances in brain function are present at RTS, which is consistent with prior studies which reported elevated resting-state functional connectivity and task-based activation after RTS (Churchill, et al., 2017b;Lovell et al., 2007). The results are surprising, as elevated scaling is thought to indicate a less taxed, more adaptive brain (Barnes et al., 2009;Churchill et al., 2016). The elevated c 1 values may therefore reflect adaptive changes that serve to maintain brain function postconcussion, similar to functional hyperconnectivity often reported after TBI (Hillary & Grafman, 2017). Alternatively, these elevations may reflect an improved disposition relative to "normal" controls, as the athletes had recently been declared fit for unrestricted sport participation. This is supported by the athletes having significantly lower symptom scores at RTS compared to their preseason baseline. Such an effect was also observed in a previous study, where recently cleared athletes scored higher on multiple psychological measures relative to matched athletic controls . Another potential contributor may be interruption in sport practice, that is, increased physiological rest. We believe it unlikely to be the primary driver of the observed effects, however. Despite substantial variability in time to RTS, the brain regions with significantly elevated c 1 at RTS did not show significant effects of CS, which is highly correlated with recovery time; nevertheless, further work is needed to validate these conclusions. Irrespective of the underlying mechanism, concussion-related abnormalities seen at RTS gradually dissipated over subsequent sessions, indicating that concussionrelated disturbances in BOLD dynamics had largely dissipated by 1 year post-RTS.
The effects of concussion on BOLD fMRI scaling values were mainly observed in posterior brain regions implicated in visual function and processing, with less frequently identified effects in temporal and thalamic regions. This has significant implications, particularly given the frequent reporting disturbances in visual function following concussion and mild TBI (Brahm et al., 2009;Capó-Aponte, Urosevich, Temme, Tarbett, & Sanghera, 2012). However, an alternative explanation for the spatial distribution of effects is that the dense vascularization and relative functional homogeneity within these brain regions may better facilitate the accurate estimation of subtle changes in scaling behavior for BOLD signals. The localization of concussion effects in this study are nevertheless consistent with a meta-analysis conducted by (Eierud et al., 2014), which showed that mild TBI is associated with enhanced posterior functional brain connectivity, supporting a distinct pattern of functional response in these brain regions after a concussion.
The analysis of CS1 and CS2 supported our third study hypothesis, as the course of brain recovery for c 1 showed significant effects of symptom severity and time to RTS. Based on the analysis of CS1, athletes with greater acute symptom burden and prolonged recovery (CS1 > 0) had elevated c 1 at RTS followed by a decline 1 year later. As the differences in brain response in this subgroup are observed at medical clearance, it remains unclear whether they represent greater adaptive neural response to injury, greater elevation in mood at RTS due to the protracted nature of recovery, or a consequence of greater physiological rest. The significant brain regions, which include inferior frontal, anterior midcingulate and rostral thalamic ROIs, are consistently implicated in motor control (Hoffstaedter et al., 2014;Lissek et al., 2007;Peyron et al., 2007). These findings suggest that motor function may be a key aspect of recovery in athletes with greater overall severity of clinical outcomes. This hypothesis is supported by a previous study of concussed athletes, which found that athletes with longer time to RTS had significantly elevated functional connectivity of motor networks at RTS (Churchill, et al., 2017b). Nonetheless, given the spatial sparsity of results, further research is needed to replicate and validate these study findings.
For the analysis of CS2, athletes with a greater acute symptom burden but relatively fast recovery (CS2 > 0) tended to have reduced c 1 at acute injury, while those with lower symptoms but a relatively slow recovery (CS2 < 0) had the opposite response. In both cases, the effects were only present at acute injury and had normalized to nearcontrol values at the time of subacute injury. The significant brain regions, which included subgenual anterior cingulate and amygdala ROIs, are both involved in emotion regulation (Etkin, Egner, & Kalisch, 2011;Haas, Omura, Constable, & Canli, 2007;LeDoux, 2003). Hence, acute disturbances in BOLD scaling for brain regions involved in emotion regulation may be an early indicator of atypical recovery time, relative to what is anticipated from symptom assessments. A possible interpretation is that athletes with CS2 > 0 have high initial worry about postconcussion outcomes, with a tendency toward greater overall symptom reporting, while those with CS2 < 0 have less worry and tend to underreport symptoms. This hypothesis is supported by prior literature showing reduced BOLD scaling in individuals with greater self-reported worry and anxiety (Churchill et al., 2015). Given the complex interplay between concussion pathophysiology, anxiety, and emotion (Mainwaring, Hutchison, Camper, & Richards, 2012;Sandel, Reynolds, Cohen, Gillie, & Kontos, 2017), future work should combine functional neuroimaging with more detailed assessments of emotional state from acute injury to RTS, in order to validate these study findings.
In terms of clinical utility, our understanding of BOLD dynamics and their relationship with concussion as a clinical syndrome remain in its early stages. Nevertheless, this study provides new insights into cerebral pathophysiology, and the relationship between brain recovery and clinical determinants of RTS. As previously noted, this study provides further convergent evidence of incomplete recovery of brain function at RTS (Churchill, et al., 2017b;Lovell et al., 2007), which may be relevant in the refinement of concussion management guidelines. This is also a promising tool for future investigation of sportrelated concussion, as analyses can be conducted on resting-state data and do not require the selection of seed regions and/or networks of interest, as in more ubiquitous measures of functional connectivity.
Even more promising, scale-free brain dynamics are, by definition, conserved across a range of timescales and potentially across different measures of brain function, as scale-free electrophysiological fluctuations are thought to underlie scale-free BOLD dynamics (Van de Ville et al., 2010). This suggests that the present findings may be readily translated into lower-cost and more portable devices, such as functional near-infrared spectroscopy (fNIRS) and electroencephalography (EEG), which are more easily implemented for patient assessment in a clinical setting.
This study provided evidence that the temporal dynamics of spontaneous BOLD fluctuations evolve over the course of concussion recovery. Nevertheless, study findings should be interpreted in the context of its limitations. There was some attrition of concussed athletes, which may lead to biased parameter estimates.
However, as noted in Section 2.1, drop out was not significantly dependent on demographics (age, sex, and concussion history) or clinical factors (symptom severity, time to RTS). Moreover, to mitigate bias and improve estimation efficiency, longitudinal analyses used LMMs and cross-sectional analyses used multiple imputation techniques. Further work is also needed to determine whether the present study findings generalize beyond university-level athletes.
We identified significant demographic effects on BOLD scaling within a relatively homogeneous athlete cohort, highlighting the importance of evaluating scaling behavior across a broader range of ages and other sport and nonathlete cohorts. In addition, although history of concussion was incorporated into analyses of both control and concussed athletes, this was based on selfreported history, which may be subject to errors and/or reporting bias. Future work should include more comprehensive documentation of concussion history and long-term follow-up to better understand the cumulative effects of multiple concussions on BOLD scaling.
In conclusion, these findings present the first examination of concussion recovery and BOLD fMRI multifractal scaling behavior, showing longitudinal effects related to recovery predominantly in brain regions involved in visual function and processing. These effects were most extensive at RTS, providing further evidence of ongoing functional brain changes beyond medical clearance. Future research should examine other measures of brain physiology in the concussion context, such as fNIRS and EEG, thereby assessing scaling of brain activity over different timescales. In addition, the effects of concussion on scaling should also be examined during cognitive tasks, to better understand how concussion affects the brain's ability to regulate mental effort. Overall, these results underscore the promise of such analysis techniques in future research and suggest that they should complement standard connectivity-based analyses of resting-state fMRI data.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.