The Confound of Hemodynamic Response Function Variability in Resting-State Functional MRI Studies

1. Introduction

Functional magnetic resonance imaging (fMRI) has contributed significantly to the advancement of neuroscience, psychiatry, and neurology over the past three decades [1,2,3]. While neural activity can be directly measured in vivo through invasive procedures, blood oxygenation level-dependent (BOLD) fMRI is a complex, indirect measure of neural activity [4] (Fig.1a), measuring local blood oxygenation variations in response to active neurons. Dilation and constriction of blood vessels modulate this process, which, in turn, is modulated through numerous non-neural and neural factors that are difficult to delineate [5]. The combination of factors that lie between neural activity and BOLD is the hemodynamic response function (HRF) [6,7]. HRF shape is characterized by its amplitude (response height, RH), latency (time-to-peak, TTP) and width (full-width at half max, FWHM) (Fig.1b). Representing neurovascular coupling in the BOLD signal, the HRF is modulated by several non-neural factors [5,8] such as hematocrit, variable size/density of vasculature, global magnetic susceptibilities, alcohol/caffeine/lipid ingestion, pulse/respiration differences, and partial volume imaging of veins [8,9,10,11,12,13,14,15]. The HRF shape varies across the brain and individuals [9,10,16].

The current ‘perspective’ article focuses on this HRF variability (HRFv) and its impact on fMRI data processing and subsequent outcome measures within the scope of human fMRI research. This topic is important because thousands of human fMRI studies are published each year, but a significant portion of those do not account for HRFv. Resting-state fMRI (rs-fMRI) and connectivity studies dominate that list. Thus, there is a need for more dialogue within the community on HRFv. We argue that HRFv causes a measurable impact on the BOLD time series, which, if ignored, will confound fMRI outcomes such as connectivity. We substantiate our argument with prior human HRFv research, including those focusing on rs-fMRI connectivity. We also present limited new data (to substantiate our 'perspective’) on HRFv across two important demographic variables: age and sex. Lastly, to cover structures of the central nervous system beyond the brain, we also discuss the implications of HRFv for spinal cord fMRI.

Figure 1. (a) FMRI is an indirect measure of neural activity. What stands between them is the hemodynamic response function (HRF). (b) The HRF is the BOLD response to a neural impulse. It has a peak (response height) occurring sometime after the neural impulse (time-to-peak) (HRF width = FWHM). The HRF always has this shape (biological), but these parameters are variable.

Figure 1. (a) FMRI is an indirect measure of neural activity. What stands between them is the hemodynamic response function (HRF). (b) The HRF is the BOLD response to a neural impulse. It has a peak (response height) occurring sometime after the neural impulse (time-to-peak) (HRF width = FWHM). The HRF always has this shape (biological), but these parameters are variable.

2. The problem of HRF variability

HRFv was first demonstrated in 1998 [10] and further examined in later years [9,17,18,19]. The HRF of a given brain region was identified to be different across individuals, and within a given individual, it was different across brain regions. Non-neural factors affecting the HRF is not merely the concern because a variable known HRF could be accounted for. Instead, the concern is that we do not understand these variable factors (and thus the HRF). The issue is alleviated in task fMRI studies by often modeling the HRF confound as time/dispersion derivatives of the canonical HRF in a general linear model (GLM) framework, which works well since BOLD is time-locked to an external stimulus. However, the outlook is different in rs-fMRI studies, which mostly ignore HRFv. The notion of HRF not being variable enough to be a serious confound was challenging to test, if not impossible, a few years ago because of the inability to estimate the HRF from rs-fMRI data. Necessary technical advancements (e.g., point-process theory [20,21]) have now resulted in rs-fMRI deconvolution techniques that are capable of HRF estimation from resting-state fMRI data [22,23,24,25]. We take a closer look at these techniques next.

3. HRF estimation from resting-state fMRI data

HRF estimation can be straightforward with two known quantities (fMRI and neural activity) and one unknown (HRF). This is possible with simultaneous fMRI and invasive recordings (e.g., [26,27]), which is hardly feasible in humans. Obtaining simultaneous rs-EEG/fMRI data and considering EEG as the neural input to deconvolve fMRI (using AFNI’s 3dDeconvolve [28]) is problematic because generative fMRI models do not consider scalp EEG as properly representing BOLD-inducing neural activity [4]. A viable alternative is estimating HRF latency from a hypercapnic challenge because breath-hold causes vasodilation and modulates cerebral blood flow (CBF) [29,30,31] in the absence of neural activity, allowing us to measure vascular latency. Chang et al. [32] utilized this to correct for vascular latency prior to connectivity analysis. The disadvantage is that an additional breath-hold scan is not always feasible or available. Breath-hold is prone to subjective performance and can be challenging in those with some neurological diseases [33,34]. Moreover, it only measures one aspect of HRF shape (TTP), while the entire HRF impacts BOLD.

To circumvent these concerns, an alternative is to perform blind deconvolution; that is, solve the mathematically ill-posed inverse problem of having two unknowns (HRF and neural activity) and one known (fMRI). This is feasible because the natural limits and biophysics of HRF and fMRI are well understood. HRF estimation then becomes a constrained optimization exercise. Such deconvolution techniques primarily focusing on task fMRI data [35,36,37,38,39,40,41,42] are generally robust to HRF misspecification within a narrow physiological range, but are not viable for estimating voxel-specific HRFs in the entire brain, especially in rs-fMRI data. Whole-brain HRF estimation is preferable with rs-fMRI data because, with task fMRI, a given task does not activate the entire brain uniformly [43], and even among activated voxels, the BOLD response is mostly non-uniform [44], sometimes leading to biologically implausible HRF estimates [43]. Some techniques that might, in principle, be viable for rs-fMRI have never been tested using rs-fMRI data [45,46].

This leaves us with four rs-fMRI deconvolution techniques that have been more widely adopted: (i) Wu et al.’s [22] data-driven rsHRF method [47], (ii) parametric generative state-space models proposed within the stochastic dynamic causal modeling (DCM) framework [48,49,23], (iii) physiologically informed DCM [24] (developed for task fMRI but can be extended to rs-fMRI), and (iv) Total Activation [25] (estimates an “activity-inducing signal” from BOLD, from which it is possible in principle to estimate the HRF by Wiener deconvolution), which has been applied in the brain [50,51] as well as spinal cord [52]. Notably, a substantial number of studies examining HRFv (described later) utilized Wu et al.’s technique [22], hence we describe it briefly here. The Wu et al. technique models rs-fMRI data as event-related time series, with events modeled as point-processes [53]. Then it estimates the best-fit HRF in a least-squares sense, using the BOLD time series at identified events. Finally, the latent neural time series is estimated using Wiener deconvolution from the measured BOLD and the estimated HRF [54]. This technique has been validated using simulations, non-invasive and invasive data [20,22,55,47], and has been applied in many recent papers [55,56,57,58,59,60,61,62]. HRFs separated by four weeks demonstrated moderate test-retest reliability (ICC=0.51) [63], which is impressive by current neuroimaging standards [64]. We next present the impact of HRFv on rs-fMRI data.

4. The confound of HRF variability on connectivity estimates

We focus here on functional connectivity (FC) [65] because most rs-fMRI studies investigate FC, or metrics derived from FC such as dynamic connectivity and graph measures. It is, however, notable that the HRF confound has also been investigated for effective connectivity models such as DCM [48], Granger causality [66,67,68,69,45], and multivariate dynamical system models [70,71,72]. Effective connectivity estimated from fMRI data is viable [66,73] as well as accurate [23,26,74,75] only after deconvolution.

FC studies have, however, largely ignored HRFv either with the assumption that the HRF is similar enough among brain regions and individuals or due to the unavailability of HRF estimation methods. Recent evidence suggests that these assumptions must be re-evaluated. Although researchers have been aware, since the early days of fMRI, that BOLD measures blood oxygenation and not neural activity, the magnitude of HRF confound on FC is being investigated only recently with the availability of deconvolution techniques. We take a closer look at these. It has been demonstrated that the HRF is variable [9,10,55,76,77] and recent reports suggest that ignoring it can introduce confounds in FC estimates [55,57,59,78,79]. In fact, our study reported significant HRFv in the brain’s default mode network (DMN) that confounded FC by about 14.7% [55]. This error in FC estimates due to HRFv (FC-error) was smaller for within-lobe FC (12.6%) vs. between-lobes (15.6%), perhaps due to more variable vasculature in the latter case.

The HRF is not only different across individuals, but also different across clinical populations due to impairments in various factors contributing to HRFv. Our prior work has highlighted this in autism [78], post-traumatic stress disorder (PTSD) [57], obsessive-compulsive disorder (OCD) [63], bipolar disorder and schizophrenia [80]. They also demonstrated that HRF impairments are significant enough to confound FC group differences; and sometimes the confound was of a similar order of magnitude as FC impairments in these diseases. Other labs have made similar observations. For example, another study found that HRF alterations confound FC [81]. HRFv also confounded rat FC [82]. Using variable HRFs that are person-specific (vs. fixed HRF) improved connectivity estimates [83]. Taken together, emerging evidence indicates that HRFv is concerning for rs-fMRI connectivity and clinical applications.

5. The importance of studying HRF variability

Non-invasive measurements are not as ‘clean’ as invasive ones. There is always an effort to make fMRI data as ‘clean’ as possible by maximizing relative variance from neural sources. Examples of such efforts include minimizing physiological/thermal noise through ultra-high-field imaging (7T) [84], improved acquisition sequences [85,86], and better denoising [21,87]. Today we are in a reproducibility crisis in functional imaging [88,89,90], indicating that further advancements are needed to make this technology clinically more useful. We still do not understand the substantial intra- and inter-subject variability of fMRI outcomes. Further characterization of this variability is timely to maximize the percentage variance in BOLD explained by the underlying neural activity.

Effective modeling of HRFv could potentially contribute to this effort. Like head motion or physiological noise, HRFv is an undesirable confound reducing fMRI data fidelity. We predict that minimizing HRFv in fMRI data will improve data quality and enhance clinical discovery. Concerns about the HRF confound also exist among the broader neuroscience community. Examples include the viewpoint of cellular neuroscientists [29], special issue articles on HRFv [91], studies on non-neural factors and BOLD [92], and investigating the link between fMRI and neural activity [93].

6. Clinical research and HRF variability

Evidence suggests that HRFv is relevant for clinical and geriatric research. HRFs in older adults are different from their younger counterparts. Older adults have longer TTP and shorter RH, largely due to vascular factors [94]. Such change is associated with Alzheimer’s disease as well [95]. Longer TTP is also linked to reduced intelligence [96]. Aberrant HRFs have been observed in stress [97], mild traumatic brain injury [98], aging [94,99,100], isolated cervical dystonia [101], and levels of consciousness [102,103]. Such HRF changes are concerning because they are at least partly driven by non-neural factors and can confound FC group differences. Although studying non-neural factors in brain disorders is a valid enterprise, attributing FC group differences entirely to neural activity is problematic. HRFv in neurological disorders is yet to be investigated; however, our prior work in psychiatric conditions (autism [78], PTSD [57], schizophrenia and bipolar [80]) demonstrated that HRF impairments in these conditions invariably confounds FC group differences, and the confound can sometimes be of the same order of magnitude as group differences. Taken together, HRFv has implications across a spectrum of cognitive, psychiatric and geriatric domains, and perhaps also in other cases in which HRFv has yet to be investigated.

7. Demographic variables and HRF variability

Non-neural factors that affect the HRF also differ across two highly relevant variables – age and sex. Aging causes vascular degradation; blood vessels of older adults are stiffer and less pulsatile, typically resulting in weaker BOLD responses to the same magnitude of neural activity (i.e., shorter RH) as well as longer time for peak BOLD activity (longer TTP) [94,99,100]. There is evidence for altered HRF in older adults [94,99], although the confound of HRFv on young vs. old FC group differences has not yet been studied. There is motivation to hypothesize that accounting for HRFv is essential for minimizing these vascular and other non-neural age-related confounds in rs-fMRI data. In fact, a recent report [99] noted that “vascular confounds in fMRI studies are common. Despite over 10,000 BOLD-fMRI papers on aging, fewer than 20 have applied techniques to correct for vascular effects.”

Factors affecting the HRF also differ between women and men. Men have lower CBF independent of neural activity [104,105], which affects the HRF [6,106]. Other factors that exhibit sex differences include vascular physiology [107], capillary microcirculation [108], and overall cerebral hemodynamics [109]. Ignoring these factors by assuming a fixed canonical HRF is a potential confound in rs-fMRI studies that report data from both sexes, and/or perform women vs. men group comparisons. Despite this, HRF sex differences have never been directly studied. Hence, we next present limited data in this context.

8. Results on HRF differences between sexes and their impact on connectivity

We utilized 7T rs-fMRI data from our earlier HRF study (N=47, 22F/25M, healthy adults) [55] (data made public [110]) (please refer to these publications for data details). Kindly note that since this is a perspective article, we have not presented comprehensive results, but we hope that the results herein will encourage extensive follow-up studies. Briefly, upon standard pre-processing, we extracted mean region-of-interest (ROI) time series from DMN regions defined by the Power-Petersen atlas [1] and estimated FC between all DMN ROI pairs using Pearson’s correlation [64,111]. This procedure was repeated for two separate pipelines: data with deconvolution (DC) and no deconvolution (NDC), which differed only in HRFv. We also obtained HRF parameters for each ROI during deconvolution.

We compared ROI-level HRF parameters between men and women (p<0.05, FDR corrected). T-test was used for RH, and Wilcoxon rank-sum test was used for discrete variables (TTP, FWHM). There were no RH differences, but men had significantly longer TTP and/or FWHM in four regions (Table 1). Longer TTP/FWHM in men could be due to CBF and vascular differences between sexes [104,105,106,107,108,109]. The mean values averaged across all ROIs and subjects were as follows: RH (women=5.01, men=4.77), TTP (women=5.61s, men=5.73s), FWHM (women=5.62s, men=5.79s).

Next, we computed the error in FC group difference T-stats resulting from HRF differences. For this, women vs. men FC stats were performed separately for DC and NDC pipelines, and the percentage difference in resulting statistical maps was computed (Fig.2). HRF differences between sexes resulted in a 15.4% median FC error. After thresholding (p<0.05, FDR corrected), all three identified FCs were false positives (FC of temporal_inf_L/frontal_sup_L/fusiform_L with temporal_mid_L). Although only four regions and three connections exhibited significant HRF and FC differences, respectively, it must not be inferred that the rest of the regions or connections were unaffected by HRFv. These HRF/FC differences were large enough to be detected for our sample size and data quality, but that does not mean that the difference can be ignored in the other connections. The testament to this fact is that HRF differences between sexes resulted in 15.4% median FC error with contribution even from a vast number of regions/FCs outside of the significant ones (Fig.2). Our prior work has demonstrated the same [55,57,78,80]. Taken together, a measurable portion of sex differences in DMN were attributable to HRFv.

As a secondary result not related to sex differences, we also report whole-brain-level average metrics of variability in HRF parameters within and between healthy young adult subjects, because these were not reported in our earlier publication [55], but they support our conclusions and encourage future work. Of course, for the nature of this article, these results could be more complex and extensive. Using the same FC data, we (i) computed within-subjects variability of each HRF parameter as the mean percentage difference between all ROI pairs; and (ii) computed between-subjects variability of each HRF parameter as the mean percentage difference between all pairs of subjects. We used a t-test for RH and a rank-sum test for TTP/FWHM (p<0.05, FDR corrected). Within-subjects, we found on average 14.1% RH variability, 13.5% TTP and 13.4% FWHM variability. L/R orbitofrontal had the largest HRF difference with other regions (86%/65% RH, 43%/24% TTP, 42%/24% FWHM variability), perhaps due to susceptibility [11] and/or vascular [99] differences. Temporal lobe had larger variability than other lobes in RH (p=0.034), TTP (p=0.039) and FWHM (p=0.037). Between-lobe variability was higher than within-lobe in RH (p=0.039), TTP (p=0.028) and FWHM (p=0.041) (possibly due to vastly different vasculature across lobes). Between-subjects HRFv (RH: 29.8%, TTP: 28.5%, FWHM: 28.1%) was larger than within-subjects. The occipital lobe varied significantly less across subjects than other lobes in RH (p=0.048) and TTP (p=0.045).

Figure 2. A flowchart illustrating our analysis and the final result obtained. The histogram shows the distribution of % FC error arising from men vs. women HRF differences. The median FC error was 15.4%, with a 95% confidence interval of [1.5% – 131.6%].

Figure 2. A flowchart illustrating our analysis and the final result obtained. The histogram shows the distribution of % FC error arising from men vs. women HRF differences. The median FC error was 15.4%, with a 95% confidence interval of [1.5% – 131.6%].

9. The HRF in the spinal cord

The central nervous system (CNS), which includes the brain and spinal cord, is a single continuous entity. But prior HRF literature is exclusively focused on the brain. Cord rs-fMRI studies have so far not accounted for HRFv, except for a recent study [52]. If systematic brain HRF changes in pathological conditions translate to the cord, it is concerning because the cord is clinically relevant for several neurological diseases (e.g., multiple sclerosis [112], chronic pain [113], amyotrophic lateral sclerosis [114], transverse myelitis [115], ataxia [116], and spinal cord injury [117]). Cord impairments are being discovered in other pathologies (Alzheimer’s disease [118] and cerebral palsy [119]), suggesting that more disorders could involve the cord than we currently understand. HRFv could confound cord FC impairments in these diseases. Hence, characterizing HRFv in the spinal cord is clinically relevant and novel.

10. Discussion and conclusions

Herein we described prior evidence for HRFv and its confound on rs-fMRI FC and elaborated on this research's importance and clinical relevance. Unexplained variability in BOLD is a more significant concern today than before because connectivity is used in sophisticated contexts such as dynamics [120], single-subject-level prediction [121], laminar fMRI [122], and precision medicine [123]. HRFv matters to a larger extent for all the desired ‘precision’ and fidelity expected of rs-fMRI today. With fast fMRI acquisition becoming prevalent [86,124], accounting for HRFv is even more critical [16] to determine the neural/vascular origin of fMRI timing differences. Thus far, FC error arising only from spatial HRFv has been quantified, and only in parts of the brain and in small samples [55]. We provided limited data for FC error between sexes. Further research is required to quantify HRFv and FC error across various within- and between-subject comparison scenarios and demographic variables.

In summary, measurements of FC often involve unexplained variance between 40% (ML prediction) and 70% (behavioral association) [121,123]. While the underlying prediction/association models may be statistically significant and perform above chance, there is still a sizeable unexplained variance. We, therefore, argue that the HRF confound, typically in the range of 10–30%, may explain a part of this variance and thus should be commonly accounted for during rs-fMRI data pre-processing.