Persistent post‐COVID headache is associated with suppression of scale‐free functional brain dynamics in non‐hospitalized individuals

Abstract Introduction Post‐acute coronavirus disease 2019 (COVID‐19) syndrome (PACS) is a growing concern, with headache being a particularly debilitating symptom with high prevalence. The long‐term effects of COVID‐19 and post‐COVID headache on brain function remain poorly understood, particularly among non‐hospitalized individuals. This study focused on the power‐law scaling behavior of functional brain dynamics, indexed by the Hurst exponent (H). This measure is suppressed during physiological and psychological distress and was thus hypothesized to be reduced in individuals with post‐COVID syndrome, with greatest reductions among those with persistent headache. Methods Resting‐state blood oxygenation level‐dependent (BOLD) functional magnetic resonance imaging data were collected for 57 individuals who had COVID‐19 (32 with no headache, 14 with ongoing headache, 11 recovered) and 17 controls who had cold and flu‐like symptoms but tested negative for COVID‐19. Individuals were assessed an average of 4–5 months after COVID testing, in a cross‐sectional, observational study design. Results No significant differences in H values were found between non‐headache COVID‐19 and control groups., while those with ongoing headache had significantly reduced H values, and those who had recovered from headache had elevated H values, relative to non‐headache groups. Effects were greatest in temporal, sensorimotor, and insular brain regions. Reduced H in these regions was also associated with decreased BOLD activity and local functional connectivity. Conclusions These findings provide new insights into the neurophysiological mechanisms that underlie persistent post‐COVID headache, with reduced BOLD scaling as a potential biomarker that is specific to this debilitating condition.


INTRODUCTION
The coronavirus disease 2019 (COVID-19) pandemic has created an ongoing public health crisis, with the disease potentially having longlasting health effects for infected individuals.A particular growing concern is post-acute COVID-19 syndrome (PACS), typically defined as COVID-related symptoms and neurological disturbances that persist over 12 weeks post infection (Shah et al., 2021;World Health Organization, 2021).This disorder is highly prevalent among survivors of COVID-19, with persistent symptoms reported in 30% or more of individuals, across a range of cohorts (Nalbandian et al., 2021).Although individuals with PACS report a diverse array of symptoms, recent studies have singled out headache as particularly noteworthy, given its high prevalence (Fernández-de-las-Peñas et al., 2021) and its debilitating impact on quality of life (Martelletti et al., 2021;Tana et al., 2022).It has been hypothesized that post-COVID headache is triggered by inflammation-mediated activation of the trigeminovascular system among individuals with pre-existing headache and/or a genetic predisposition to migraine (Straburzyński et al., 2023;Tana et al., 2022), although COVID-related neural injury may also play a role (Galea et al., 2022).At present, we have an incomplete understanding of the etiol-ogy of headache in PACS and its effects on neural function, making this a critical area of ongoing research.
Although understudied in the context of COVID-19, there is an extensive body of neuroimaging literature examining headache in other conditions, such as migraine, cluster headache, tension-type headache, and medication overuse headache (e.g., see Chen et al., 2017;Ferraro et al., 2012;Tu et al., 2020;Wang et al., 2014).These studies have frequently reported altered brain function within subcortical and cortical regions implicated in pain processing and modulation (Chong et al., 2019;May, 2006), including the amygdala, thalamus, basal ganglia, anterior insula, anterior cingulate, sensorimotor cortex, prefrontal cortex, and posterior parietal cortex.Much of the literature also shows changes in distributed functional networks that are not specific to headache, such as the default-mode, salience, and executive control networks (Chong et al., 2019;Maleki & Gollub, 2016), emphasizing that this condition tends to have a global impact on brain functioning.By contrast, previous functional neuroimaging studies of COVID-19 have identified alterations in brain activity associated with persistent symptoms, including reduced neurometabolism in subcortical regions (Guedj et al., 2021;Sollini et al., 2021) and patterns of altered inter-regional communication related to olfactory dysfunction (Esposito et al., 2022) and to neuropsychiatric symptoms (Cattarinussi et al., 2022).At present, however, the effects of persistent post-COVID headache on brain function remain understudied.One promising approach for assessing these effects is through the dynamics of resting brain activity.
In a healthy brain, resting neural activity exhibits scale-free behavior, in which fluctuations occur across a wide range of timescales, with no scale having a dominant role.The fluctuations are thus "self-similar" so that if a segment of a time series is taken and magnified, it statistically resembles the whole.Stated more precisely, a time series of brain activity x(t) exhibits self-similarity if, for any positive scalar value a, it is identically distributed with a dilated and scaled version of itself x(t) ∼ a −H x(at).Scaling behavior is then fully described by the Hurst exponent H, with higher values indicating more scale-free signal.In signal analysis, scaling is usually defined via a power-law relationship between frequency f and the corresponding power spectral density (PSD), expressed as PSD(f) ∼ |f| − .The scaling exponent  can then be mapped to the Hurst exponent via the relationship H = ( + 1)∕2 under a fractional Gaussian noise model (Eke et al., 2000).Other scaling estimators have also been developed, including time-domain methods (DFA; Peng et al., 1994Peng et al., , 1995) ) and wavelet-based techniques (Abry & Veitch, 1998;Abry et al., 1995).In addition, more advanced techniques have been used to describe the spectrum of time-varying scaling exponents h instead of a single exponent H (Lashermes et al., 2005;Wendt & Abry, 2007;Wendt et al., 2007).These techniques have been applied to functional neuroimaging data obtained from both healthy and patient cohorts (Ciuciu et al., 2012;He, 2011He, , 2014;;Van de Ville et al., 2010;Wink et al., 2008).Studies have found that scaling behavior is often suppressed (i.e., the exponent H is decreased) under conditions of increasing physiological and psychological demand, for example, during challenging and/or unfamiliar tasks (Barnes et al., 2009;Churchill et al., 2016;Ciuciu et al., 2012), with increased age (Churchill et al., 2016), under psychological distress (Churchill et al., 2015;Tolkunov et al., 2010), and following mild traumatic brain injury (Churchill et al., 2020).Such findings suggest that H may similarly be sensitive to PACS, in terms of the mental burden of post-COVID symptoms (Peghin et al., 2021) and any underlying neural injury.This is particularly relevant for persistent headache, which may cause significant psychological distress (Kristoffersen et al., 2018) and tends to interfere with normal cognitive function (Moore et al., 2013;Vuralli et al., 2018).
The present cross-sectional observational study examined this question using resting-state blood-oxygenation-level dependent func-

Scaling analysis
The overall goal of the analyses was to summarize the dynamics of the spontaneous, arrhythmic fluctuations of resting-state BOLD time series x(t).One approach involves analyzing signal power at different frequencies f using a PSD estimator.However, both the BOLD signal and the underlying neural activity have a broad-band PSD, in which power declines smoothly with increasing frequency (Fox et al., 2007), hence, no specific set of frequencies or time scales can be singled out for analysis.For this reason, scaling analysis is an appealing alternative: instead of examining the absolute BOLD power at a specific frequency band, this approach characterizes the relationships between relative BOLD power at different frequences in terms of power-law scaling behavior.By convention, x(t) is deemed scale invariant if the expression PSD(f) = C|f| − holds for a wide range of frequencies, where  > 0 and C is an arbitrary scaling constant.
This has non-trivial implications, as the relative spectral content at a given frequency f is then entirely dictated by a single scaling exponent .This parameter is also related to the Hurst exponent by H = ( + 1)∕2 for a fractional Gaussian noise model (Eke et al., 2000), where H is considered an index of temporal dependence.A value of 0.5 < H < 1.0 indicates positive long-range autocorrelation in the signal of interest, where a high (or low) value at x(t) tends to be followed by high (or low) values at x(t + 1) and subsequent x(t + n).Signals with an H value near 1.0 have highly persistent autocorrelations for time lags of n ≫ 1.
Conversely, signals with an H value near 0.5 show autocorrelations that decay rapidly for n ≥ 1.Values of 0 < H < 0.5 may also be observed and denote persistent long-range anticorrelations, where a high (or low) value at x(t) tends to be followed by low (or high) values at x(t + 1), with switching behavior extending over time lags of n ≫ 1.In practice, PSDbased scaling analysis involves plotting log(PSD(f)) values against log(f) values and calculating the linear slope coefficient  from the resulting scatter plot, which is then used to estimate H.
For the present study, the PSD was calculated at the voxel level using a Welch estimator (Hanning windowed, with 50% overlap between windows), with band integration over 0.02 Hz intervals to improve estimation stability.The scaling coefficient was then obtained by plotting log(PSD(f)) against log(f) in the range of 0.015 to 0.225 Hz, followed by ordinary least squares estimation of the slope coefficient .After verifying the appropriateness of the fractional Gaussian noise model (Eke et al., 2000) and linear goodness of fit, the Hurst coefficient H was obtained for each voxel.Further testing also affirmed that the scaling relationships identified in this study were not substantially affected by a change in the frequency lower bound (tested: 0.015, 0.02, 0.025, 0.03 Hz), nor the upper bound (tested: 0.225, 0.200, 0.1875, 1750 Hz).
While the main analyses of this study used a PSD-based estimator of scaling behavior, numerous other scaling estimators exist.Given the ongoing debate surrounding the advantages and drawbacks of different methods in the context of BOLD fMRI, we have also examined BOLD scaling behavior using three other well-established techniques that have distinct estimation procedures.They include two "monofractal" approaches that describe scaling behavior in terms of a single exponent H, along with a "multifractal" approach that summarizes the data in terms of a spectrum of time-varying scaling exponents.The techniques are detrended fluctuations analysis (DFA; Peng et al., 1994Peng et al., , 1995)), wavelet monofractal analysis (WMA; Abry & Veitch, 1998;Abry et al., 1995), and wavelet leader multifractal analysis (WLM; Lashermes et al., 2005;Wendt & Abry, 2007;Wendt et al., 2007).Details of the estimation procedures and the subsequent analysis results are presented in Appendix 2 in the Supporting Information.

Analysis of clinical and demographic data
Participant demographics are listed in A series of two-sample analyses then tested for group differences in demographic variables.The mean differences between groups were calculated, along with bootstrapped standard errors, bootstrap ratios (BSRs; z-scored statistics of effect, based on the ratio of mean/SE) and two-tailed percentile p-values; notable differences were identified at an uncorrected threshold of p < .05.For the COVID-19 group, headache was also compared to other non-headache symptoms.The frequencies of other ongoing symptoms was reported, along with their bootstrapped 95%CIs, and frequencies were compared to headache with reporting of the mean differences, bootstrapped 95%CIs and twotailed p-values.Correlations were also calculated between ongoing headache and the presence of other ongoing symptoms, with reporting of Spearman correlation coefficients, bootstrapped 95%CIs and two-tailed p-values.

Global effects of COVID-19 on BOLD scaling
For each participant, the global Hurst exponent H glob was calculated as the average overall gray matter voxels.To assess for the effects

Regional effects of COVID-19 on BOLD scaling
Subsequent analyses sought to localize the regional effects of COVID-19 and post-COVID headache on H values, for the group contrasts yielding significant effects on global scaling estimator H glob .For each voxel in the brain, a GLM was fitted on local H values in the same manner as the previous section, obtaining coefficient contrasts, along with bootstrapped 95%CIs, BSRs, and two-tailed p-values.Then, for each paired group contrast, significant voxels were identified by thresholding at an uncorrected p-value threshold of .005,followed by cluster-size thresholding at p = .05,using the AFNI 3dFWHMx program to estimate spatial smoothness and the AFNI 3dClustSim program to estimate the corresponding minimum cluster-size threshold.For significant brain voxels, the regional effect sizes were then reported in terms of the BSR values and cluster summary reports were also generated.

Comparison with alternative BOLD measures
To provide further context for the observed group differences in scaling behavior, more conventional measures of resting-state BOLD signal were also examined including the amplitude of low-frequency fluctuations (ALFF), which computes total BOLD power summed over the range of 0.015-0.08Hz, and its normalized equivalent, the fractional ALFF (fALFF), which normalizes this power sum by the total spectral power.Both metrics were calculated using the same PSD method as in the scaling analyses.Note: For each symptom, prevalence is given as the percentage of the overall coronavirus disease 2019 (COVID-19) group (N = 57) that report ongoing issues, along with a bootstrapped 95% confidence interval (95%CI).The mean difference in prevalence relative to headache is also computed, with a 95%CI and percentile-based p-value.For symptoms with prevalence ≥ 10% (and thus stable bootstrap intervals), Spearman correlations with ongoing headache are computed, with bootstrapped 95%CIs and p-values.
These measures were calculated for a 27-voxel cubic ROI placed on the center of mass for each cluster identified as having significant group differences in scaling effects in the section above, with subsequent averaging over all clusters to obtain a single set of (ALFF, fALFF, Lconn, Gconn) values per participant.Afterward, GLM analyses were conducted to evaluate group differences for each of the BOLD measures as in the previous sections, obtaining coefficient contrasts, along with bootstrapped 95%CIs, BSRs, and two-tailed p-values.Although there was substantial reporting heterogeneity between individuals, the COVID-H+ group tended to report higher frequencies of cognitive issues and body pain/ache relative to other subgroups.No other noteworthy effects were found for any of the other group comparisons and/or symptoms, with all |BSR| ≤ 1.72 and p ≥ .074.

Regional effects on scale-free dynamics
Figure 2 displays regional H values and between-group differences.Note that the contrast between COVID-H-and control groups was not plotted, as no significant differences were identified.Standardized effect sizes are displayed as z-distributed bootstrap ratio (BSR) statistics.(c) Consensus areas are depicted where COVID-H+ and COVID-Hr groups show significant differences relative to both controls and the COVID-H-group.

DISCUSSION
The COVID-19 group without headache did not significantly differ from the control group in terms of global or local H values.This suggests that, contrary to study hypotheses, PACS alone does not substantially alter functional brain dynamics relative to symptomatic non-COVID infection, despite evidence of COVID-19 targeting brain tissues (Baig et al., 2020;Galea et al., 2022).In terms of the psychological burden, results are also consistent with the similar levels of symptom burden seen in this study's COVID-19 and control groups as examined in Churchill et al. (2023).However, these findings may be due to the relatively mild sequelae in the current sample of non-hospitalized individuals with COVID-19, compared to others who may experience flagrant PACS.Given the high rates of high levels of distress often reported among individuals with PACS (Leviner, 2021) Churchill et al., 2016), with greater effects seen in more challenging conditions relating to task difficulty, task unfamiliarity, and participant age (Churchill et al., 2016).Hence, the present study findings may reflect the greater demands of post-COVID headache related to, for example, sensory and pain processing.Alternatively, previous studies have identified H decreases due to psychological distress (Churchill et al., 2015) and abnormal BOLD scaling for individuals with high trait anxiety (Tolkunov et al., 2010).Given that symptoms of mental distress are common in headache disorders (Kristoffersen et al., 2018), this may also contribute to the effects seen in the present study.Finally, a study of mild traumatic brain injury found decreased H during the acute phase of injury (Churchill et al., 2020), along with chronic declines in H among the most symptomatic individuals.Similar to PACS, this cohort is likely experiencing a combination of diffuse neural injury and psychological distress relating to post-concussion symptoms.More broadly, scale-free dynamics are a property of many biological systems, with links to concepts such as complexity and criticality (Delignières & Marmelat, 2012).Reduced scaling behavior in biological systems generally signifies a loss of system complexity (Goldberger et al., 2002) and thus a more constrained, less adaptive state, which may be vulnerable to further disruption.
Unexpectedly, the COVID-19 group who had recovered from headache not only had elevated global H values relative to COVID-19 with ongoing headache but also had slightly elevated values relative to the non-headache COVID-19 and control groups.It is interesting to note that a similar finding was observed in the previously cited study of mild traumatic brain injury (Churchill et al., 2020), where individuals showed elevated H values at medical clearance.There, it was hypothesized that this represents a relative improvement in mental state, compared to uninjured controls, due to having recently recovered.Given the often debilitating nature of post-COVID headache (Tana et al., 2022), it is possible that a similar effect is seen here; future studies examining perception of post-COVID headache and recovery may be informative in interpreting these findings.
Further headache-related hypotheses posited decreased H primarily in areas previously identified across different forms of headache disorder (May, 2006), including the amygdala, anterior cingulate, anterior insula, medial prefrontal, and sensorimotor regions.The present findings again partly support this hypothesis, with alterations in insular and sensorimotor regions reliably identified.The contrast of ongoing versus recovered headache subgroups revealed further effects in the anterior cingulate and amygdala, although the prefrontal cortex was not identified.These areas have been routinely observed in studies of migraine, cluster headaches, and medication overuse headache (Chong et al., 2019).The identified regions are also routinely associated with pain and noxious stimuli, with the somatosensory cortex and amygdala being among the most frequently reported (Iannetti & Mouraux, 2010).
Overall, these results support the interpretation that reduced H values seen in the COVID-19 group with headache are at least partly driven by the increased sensory and pain processing demands.
Further analyses using other estimators of the scaling exponent H also indicated that the observed effects generalized across all of the tested estimators, including DFA, time-frequency monofractal meth-ods (WMA) and time-frequency multifractal methods (WLM).This affirms that the present findings related to post-COVID headache are not unique to the PSD methodology.This finding is also consistent with prior studies of BOLD scaling behavior, in which the main study conclusions were found to generalize across different estimators of BOLD scaling behavior (Churchill et al., 2015(Churchill et al., , 2016)).
Additional comparison of H with other measures indicated that areas of significantly reduced (or increased) BOLD signal scaling associated with post-COVID headache had similar reductions (or increases) in local BOLD activity and functional connectivity.This is consistent with prior studies showing correlations of H with indices of BOLD variability and functional connectivity strength in healthy adults (Churchill et al., 2016;He, 2011).The present results provide further evidence that a change in scaling dynamics of the brain is reflected in other aspects of brain function.For example, the ALFF and fALFF indices reflect the "dynamic range" of the neurovascular response (Yang et al., 2007;Zou et al., 2008), whereas Lconn reflects the local homogeneity of BOLD response and is related to segregative brain function (Zang et al., 2004).Interestingly, Gconn, which reflects global integration of brain regions, does not show significant decreases in H, indicating the changes in BOLD scaling during ongoing headache may be more relevant to local neural functioning in this cohort.
This study had some limitations that should be acknowledged, as they qualify the interpretation of the results.First, the fMRI acquisition had a run length of ∼9 min, which is consistent with standard resting-state protocols and balances the demands of adequate data and patient wakefulness (Birn et al., 2013;Tagliazucchi & Laufs, 2014).
However, physiological time series on the order of 10 3 or more are preferred for robust scaling analysis (Eke et al., 2000); future research employing multiple runs may be able to bypass these issues and explore scaling behavior over longer time scales.The sample sizes are also somewhat small for participant subgroups, although given the paucity of data examining PACS and post-COVID headache, this remains a valuable dataset.In particular, the control group is relatively unique, having non-COVID infection and similar post-infection symptom profiles, compared to the more common literature practice of contrasting individuals with COVID-19 against uninfected controls.
The bootstrapped regression approach was chosen to preserve power in small and unbalanced samples, while the large effects sizes, robustness across scaling measures, and predefined hypotheses provide some internal validation of study findings.Nevertheless, it will be critical for tional magnetic resonance imaging (BOLD fMRI) data collected as part of the Toronto-based NeuroCOVID-19 study(MacIntosh et al., 2021).The study compared whole-brain functional dynamics of self-isolating individuals who tested positive for severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and subsequently experienced persistent lingering symptoms, relative to controls who had cold or flulike symptoms but tested negative for SARS-CoV-2, with both groups imaged an average of 4-5 months after COVID testing.Analysis was conducted after dividing the COVID-19 group into those who had no headache symptoms (COVID-H-), those who had ongoing headache symptoms (COVID-H+), and those whose headache symptoms were resolved at the time of imaging (COVID-Hr).These groups were compared in terms of the global mean H values and regional H values, to enable testing of the hypothesis that post-COVID-19 condition produces a graded effect on H, with COVID-H-showing reductions relative to controls, COVID-H+ showing further reductions from the COVID-H-, and COVID-Hrshowing similar values to COVID-H-.The study design also enabled testing of a second, related hypothesis that the brain regions implicated in COVID-H+ include consensus regions identified in studies of headache disorder(Chong et al., 2019), i.e, the amygdala, thalamus, basal ganglia, anterior insula, anterior cingulate, sensorimotor cortex, prefrontal cortex, and posterior parietal cortex.Finally, supplemental analyses characterized the relationship between alterations in H and other commonly used indices of resting brain function, including BOLD activity in the 0.015-0.08Hz frequency band, along with indices of local and global functional connectivity.
of COVID-19 and post-COVID headache on H glob , a general linear model (GLM) was fitted on the data from all groups (control, COVID-H-, COVID-H+, COVID-Hr) with three binary regressors denoting membership in the different COVID-19 subgroups; the GLM also adjusted for effects of age and sex.Coefficients of effect were obtained for each group, along with their bootstrapped distributions.Contrasts were then obtained between all group pairs, with bootstrapped 95%CIs, BSRs, and -tailed p-values.Groups with significantly different H glob values were identified after adjusting for multiple comparisons at a false discovery rate (FDR) of 0.05.Sensitivity analyses were also con-ducted, examining the robustness of analysis results under different GLM modeling choices, including the effects of incorporating supplemental demographic and clinical variables; the results are summarized in Appendix 3 in the Supporting Information.

Figure
Figure 2a plots the mean H map, computed over all participants.Values tended to be highest in the posterior occipital, temporal, and parietal regions, along with the cingulate cortex.Intermediate values are observed frontally, and lower values are observed in subcortical and medial temporal areas.Figure 2b plots brain regions showing significant differences in local H.In addition, Figure 2c shows regions where COVID-H+ and COVID-Hr groups show significant differences relative to both controls and the COVID-H-group, with clusters summarized in Table 3.Only the comparison of COVID-H-and control
future research to validate and replicate the current results.Future work should also collect more detailed information about headache symptoms, including phenotype (e.g., migraine vs. tension type), severity, and frequency, in order to gain a more nuanced understanding of the relationship between BOLD dynamics and headache.Finally, while the present study excluded participants with prior neurological conditions, future research should examine individuals with PACS and pre-existing headache disorder.There is emerging evidence that the course of headache is significantly altered during COVID-19 infection(Waliszewska-Prosół & Budrewicz, 2021), and neuroimaging of these cohorts may provide greater insights into the mechanisms underlying post-COVID headache.

Table 1
Summary of demographic and clinical data for study participants.
normality of demographic data via skewness and kurtosis permutation testing at a threshold of p< .05(5000 resamples), means and standard deviations were reported for measures that did not deviate from normality, and medians with upper and lower quartiles were reported for those that did.All participants completed a questionnaire evaluating symptom status for nine items: fever, cough, sore throat, shortness TA B L E 1 had occurred but resolved, or (3) was currently ongoing.The present study focused on "headache" as the main symptom of interest.The frequency of ongoing and resolved headache symptoms was reported, along with bootstrapped estimates of the 95% confidence intervals (95%CIs).Note that all bootstrap statistics reported in this study are calculated based on 2000 resampling iterations.
Summary of non-headache symptom data for study participants.
Measures of BOLD functional connectivity were also obtained, including the local connectivity (Lconn), obtained by calculating the pairwise Pearson correlation between all voxel pairs within a region of interest (ROI) and taking the average value, and the global connectivity (Gconn), obtained by taking the mean BOLD time series in an ROI and calculating the mean absolute connectivity with all other voxels in the mask of cerebral grey matter.TA B L E 2