Resting-state fMRI signals contain spectral signatures of local hemodynamic response timing
Abstract
Functional magnetic resonance imaging (fMRI) has proven to be a powerful tool for noninvasively measuring human brain activity; yet, thus far, fMRI has been relatively limited in its temporal resolution. A key challenge is understanding the relationship between neural activity and the blood-oxygenation-level-dependent (BOLD) signal obtained from fMRI, generally modeled by the hemodynamic response function (HRF). The timing of the HRF varies across the brain and individuals, confounding our ability to make inferences about the timing of the underlying neural processes. Here, we show that resting-state fMRI signals contain information about HRF temporal dynamics that can be leveraged to understand and characterize variations in HRF timing across both cortical and subcortical regions. We found that the frequency spectrum of resting-state fMRI signals significantly differs between voxels with fast versus slow HRFs in human visual cortex. These spectral differences extended to subcortex as well, revealing significantly faster hemodynamic timing in the lateral geniculate nucleus of the thalamus. Ultimately, our results demonstrate that the temporal properties of the HRF impact the spectral content of resting-state fMRI signals and enable voxel-wise characterization of relative hemodynamic response timing. Furthermore, our results show that caution should be used in studies of resting-state fMRI spectral properties, because differences in fMRI frequency content can arise from purely vascular origins. This finding provides new insight into the temporal properties of fMRI signals across voxels, which is crucial for accurate fMRI analyses, and enhances the ability of fast fMRI to identify and track fast neural dynamics.
Editor's evaluation
This manuscript addresses the important issue of hemodynamic response function (HRF) variability across brain areas and will be valuable to researchers who use fMRI and other types of functional imaging that rely on neurovascular coupling. Using simulations and experiments, the authors provide compelling evidence that differences in the HRF can impact spectrum-based metrics such as ALFF and fALFF. A better understanding of the variability of the HRF is critical for the proper interpretation of activation onset times and of differences observed in clinical populations where both neural and vascular alterations can be expected.
https://doi.org/10.7554/eLife.86453.sa0eLife digest
Functional magnetic resonance imaging (fMRI) is a tool that can be used to non-invasively measure the activity of the human brain. Active parts of the brain require more oxygen, which increases blood flow to these areas. fMRI can detect these changes, and its signal reflects the coupling between brain activity and changes in blood flow.
The mechanism that couples brain activity to blood flow is known as the ‘hemodynamic response’, and its timing varies across the brain. Therefore, to interpret fMRI signals correctly and use them to measure underlying brain activity, it is necessary to understand how the response changes across the brain.
Current methods for probing hemodynamic response variation are either limited to specific brain regions or require patients to hold their breath – something not all groups of patients can do. To solve this problem, Bailes et al. investigated whether resting-state fMRI signals contain information about how the hemodynamic response changes across the brain. This information could then be used to better infer brain activity from fMRI measurements.
The experiments showed that resting-state fMRI signals can be used to characterize and predict the timing of the hemodynamic response. Specifically, the frequencies in resting-state fMRI signals are impacted by changes in the hemodynamic response and can therefore be used to predict hemodynamic timing. Additionally, Bailes et al. showed that these predictions are better than those obtained in experiments requiring patients to hold their breath, which is the current gold standard. The findings also demonstrate that the information from the frequencies of resting-state fMRI signals should be interpreted carefully, as differences in these frequencies can have a non-neural origin.
Bailes et al. propose a highly generalizable approach for mapping and predicting variations of the hemodynamic response across the whole brain. These findings provide insights into the time-related properties of fMRI signals that are crucial for accurate analyses. This will be of particular importance as the field moves towards fMRI studies focused on rapid neural dynamics and higher-level cognition.
Introduction
Functional magnetic resonance imaging (fMRI) enables non-invasive measurement of human brain activity via the hemodynamic response. When activity in a population of neurons changes, these changes give rise to the blood-oxygenation-level-dependent (BOLD) signal measured in most fMRI studies (Ogawa et al., 1990). Thus far, however, BOLD fMRI has exhibited relatively limited ability to provide the fine-grained temporal information necessary for deepening our understanding of brain dynamics. This is due to the fact that the signals obtained from BOLD fMRI are not direct measures of neural activity, but rather reflect the coupling between neuronal activity and the hemodynamic response, which evolves on a time course of seconds (Kwong et al., 1992; Ogawa et al., 1990). This coupling between neural activity and the BOLD signal can be represented by the hemodynamic response function (HRF) (Aguirre et al., 1998; Handwerker et al., 2004). The properties of the HRF depend on many interconnected factors, including the effects of local vascular architecture and cerebrovascular dynamics, that vary substantially across the brain and between individuals (Aguirre et al., 1998; Handwerker et al., 2004; Len and Neary, 2011; Logothetis et al., 2001). The relative timing and shape of the HRF, therefore, also varies considerably across brain regions and even between neighboring voxels (Buckner et al., 1998; Lee et al., 1995; Birn et al., 2001; Siero et al., 2011; Miezin et al., 2000; Siero et al., 2015). Hemodynamic response temporal lag variation is substantially larger than many neural effects of interest, introducing variability on the order of several seconds (Buckner et al., 1998; Lee et al., 1995; Birn et al., 2001; Siero et al., 2011). Thus, to enable inferences about the relative timing of neural activity using signals obtained from BOLD fMRI, it is crucial to understand the variations in HRF timing across the brain.
Advances in acquisition technology now allow high-resolution whole brain fMRI data to be acquired at fast (<500ms) rates (Polimeni and Lewis, 2021; Chen et al., 2019; Barth et al., 2016; Chiew et al., 2018; Hennig et al., 2007; Setsompop et al., 2016), suggesting that fMRI could provide a unique tool to noninvasively track temporal sequences of neural activity (Menon et al., 1998) across the entire brain. Indeed, recent studies have revealed highly structured temporal dynamics using fMRI and suggest that fMRI can enable whole-brain mapping of temporal sequences (Lee et al., 2013; Setzer et al., 2021; Raut et al., 2021; Mitra et al., 2016; Vizioli et al., 2018). Furthermore, hemodynamic signals have been shown to contain more information about fast and high-frequency activity than previously thought (Lee et al., 2013; Smith et al., 2012; Chen and Glover, 2015; Lewis et al., 2016; Sasai et al., 2021). Studies examining individual brain regions have demonstrated that fMRI can achieve impressive temporal precision within regions, on the order of 100ms (Lin et al., 2013), meaning that high fidelity temporal information is present within these hemodynamic signals. However, a key remaining challenge is that the hemodynamic differences across the brain confound our ability to infer the timing of the underlying neural activity from BOLD fMRI. Specifically, if a given brain region shows earlier BOLD activity, it could be due to faster neural activity or simply due to a faster hemodynamic response in that region. Fully exploiting the higher temporal resolution provided by fast fMRI techniques will therefore ultimately require accounting for differences in the temporal dynamics of the hemodynamic response across the whole brain.
Despite the well-known heterogeneity of hemodynamic timing, most common analysis approaches for BOLD fMRI data assume a standard, canonical HRF shape throughout the brain (Glover, 1999). This approach is understandable, since the true HRF is not known, but it nevertheless cannot account for the vascular confound introduced by hemodynamic response variability. Incorrect assumptions about the shape and timing of the HRF can lead to incorrect inferences regarding the underlying neural activity (Gonzalez-Castillo et al., 2012; Lindquist et al., 2009; Handwerker et al., 2012), and studies that assume a whole brain canonical HRF are unable to decouple the neural and vascular components of the BOLD signal (Handwerker et al., 2004; Rangaprakash et al., 2018; Rangaprakash et al., 2017; Deshpande et al., 2010; Chang et al., 2008). Even when using flexible modeling approaches, such as basis sets or finite impulse response models, it is not possible to determine whether a given region’s faster fMRI response reflects fast neural activity, or simply faster local neurovascular coupling (Handwerker et al., 2012).
The fact that most studies do not account for variations in HRF dynamics is largely due to methodological challenges. Previous work has demonstrated that it is possible to quantify hemodynamic lags across brain regions, and even on a voxel-wise level, to detect the relative order of BOLD responses with high temporal precision (Siero et al., 2011; Miezin et al., 2000; Siero et al., 2015; Lin et al., 2013; Chang et al., 2008; Sicard and Duong, 2005; Lin et al., 2018; Misaki et al., 2013; Kastrup et al., 1999; Liu et al., 2017; Bright et al., 2009; Chen et al., 2021; Wu et al., 2013; Wu et al., 2021). One such method is the use of a stimulus paradigm that drives activity in particular brain regions where the neuronal response properties are relatively well understood and controlled such as primary sensory or motor cortices (Handwerker et al., 2004; Birn et al., 2001; Siero et al., 2011; Miezin et al., 2000; Lin et al., 2018). However, this approach cannot be applied to the majority of the brain, where the neuronal response properties are not known ahead of time. Alternatively, using a breath hold or similar hypercapnic challenge can modulate cerebral blood flow (CBF) to all vascularized regions with minimal changes in cerebral metabolic rate of oxygen (CMRO2), allowing mapping of vascular latencies (Chang et al., 2008; Sicard and Duong, 2005; Kastrup et al., 1999; Liu et al., 2017; Bright et al., 2009; Chen et al., 2021; Pinto et al., 2021). However, breath hold tasks are not suitable for all subject populations, as some patients may have difficulty complying with the breath hold task. Furthermore, breath hold tasks modulate and measure cerebrovascular reactivity (CVR), which contributes to neurovascular coupling but is a distinct process. While neurovascular coupling reflects the alterations in local hemodynamics that occur in response to changes in neural activity, CVR is specifically a measure of a blood vessel’s capacity to dilate and constrict in response to a vasoactive stimulus, and does not include the extensive metabolic and molecular factors that also drive neurovascular coupling (Pinto et al., 2021; D’Esposito et al., 2003; Iadecola, 2017). In fact, there is evidence that while CVR is affected by healthy aging, some metrics of neurovascular coupling are not, hinting that distinct mechanisms may shape these two patterns (Stefanidis et al., 2019; West et al., 2019; Grinband et al., 2017).
A task-free, neurovascular-based approach for detecting the lags of intrinsic neurovascular coupling would therefore be broadly relevant for analyzing fMRI data. A potential alternative route towards identifying local hemodynamic properties is to examine the properties of resting-state fMRI data. Resting-state fMRI signals reflect neurovascular coupling induced by spontaneous neural activity (Mateo et al., 2017; Ma et al., 2016) and confer the additional benefit of being task-independent, which makes it a viable scan type for patient populations. Moreover, unlike stimulus- or task-based paradigms for mapping local hemodynamic response timings in specific brain regions, resting-state fMRI can be used to examine the HRF across the whole brain. These advantages have prompted past research into the utility of resting-state fMRI signals to estimate the HRF itself with prior work utilizing deconvolution approaches to explore HRF timing in resting-state data (Wu et al., 2013; Wu et al., 2021; Sreenivasan et al., 2015). However, these approaches require assumptions about the underlying neural events, which are not known. We therefore investigated whether intrinsic signatures of local neurovascular coupling dynamics are present in the resting-state fMRI signal.
Our goal was to understand whether information about the temporal dynamics of the hemodynamic response is present in resting-state fMRI data. We first used simulations of the BOLD response to illustrate how distinct, physiologically relevant HRF shapes should produce marked differences in the frequency content of resting-state signals. Next, we verified this result in fast fMRI data collected at 7 Tesla (T) using visual stimulation to induce a neural response with known timing in primary visual cortex (V1). We quantified the temporal delay of voxels in V1 in response to a controlled, oscillating visual stimulus, and found that voxels with fast and slow hemodynamic responses exhibited distinct resting-state spectral features. We further extended our analyses to the visual thalamus (lateral geniculate nucleus, LGN) and found that this principle generalized to subcortex. To understand the potential of this information as a tool to predict the temporal properties of individual voxels, we then trained classifiers to use information from the resting-state spectrum to classify voxels as being fast or slow cortical voxels, or even faster LGN voxels. We found that resting-state signals were better predictors of voxel-wise differences in relative hemodynamic timing than latencies measured from a gold standard breath hold task. Our results establish that information about hemodynamic timing can be extracted from the frequency spectrum of resting-state fMRI signals. This demonstrates that resting-state fMRI can provide a way to understand and predict the temporal dynamics of the HRF across the brain, which is critical for interpreting neural activity using BOLD fMRI.
Results
The temporal dynamics of the HRF profoundly impact the spectrum of simulated BOLD responses
Previous modeling work has illustrated that narrower HRFs should result in BOLD responses containing more high frequency power (Chen and Glover, 2015), suggesting that local variations in HRF timing should manifest as local variations in the frequency content of BOLD signals. We therefore hypothesized that the frequency spectrum of BOLD dynamics in the resting-state can be used to infer the relative timing of the task-driven hemodynamic response. We first aimed to illustrate this property by simulating the BOLD frequency response using HRFs with different temporal dynamics. If we assume that the BOLD response is a linear time-invariant system, we can compute the BOLD response as a convolution between a given input (i.e. the stimulus) and the characteristic input response of the system (i.e. the HRF). Then, by varying the frequency of the stimulus we can construct a spectrum of the BOLD frequency response for HRFs with faster or slower dynamics (Figure 1A). We performed this simulation using six different HRFs (Figure 1B) with physiologically representative values for their time-to-peak (TTP), full width at half maximum (FWHM), and amplitude (Siero et al., 2011). We observed that while HRFs with faster dynamics produced less power in the low frequency bands compared to those with slower dynamics, they also showed a shallower decline in power at higher frequencies (Figure 1C). Furthermore, this effect was preserved when we normalized the different HRFs to have the same peak amplitude (Figure 1—figure supplement 1A), demonstrating that this phenomenon is not solely due to the higher amplitude of slower HRFs (Figure 1—figure supplement 1B). This effect was also preserved if we accounted for the 1/f-like spectral pattern that neural activity displays (Figure 1—figure supplement 1C–E). And finally, although prior work has shown that the TTP and FWHM of the HRF are correlated (Siero et al., 2011; de Zwart et al., 2005), we also tested the effect of holding one parameter constant (Figure 1—figure supplement 2). Varying the FWHM had a more profound impact on the frequency spectrum compared to varying the TTP (Figure 1—figure supplement 2), as expected due to the higher frequency content of narrower HRF shapes, although short TTPs also had a smaller effect on the spectrum. These simulations demonstrate that there should be profound differences in the relative power at low versus high frequencies for voxels with fast vs. slow HRFs, which we hypothesized can be quantified and related to the temporal dynamics of the hemodynamic response.
Features of the resting-state spectrum show significant differences between fast and slow voxels
Based on these simulation results, we then examined the frequency content of resting-state fMRI data. Spontaneous BOLD fluctuations captured in resting-state fMRI are linked to fluctuations in neuronal activity (Mateo et al., 2017; Ma et al., 2016; Fox and Raichle, 2007), and accordingly, reflect the neurovascular coupling mechanisms that link neural fluctuations to BOLD fluctuations. To empirically test the prediction that fast and slow voxels should have distinct frequency content in the resting-state, we first used a task paradigm to identify voxels in V1 with consistently fast or slow hemodynamic responses. To drive continuous oscillations in V1 we presented the subjects with a 12 Hz counterphase flickering radial checkerboard with the luminance contrast of the checkerboard modulated in time as a sine wave of 0.05 Hz (Figure 2A, top). We used a combination of an anatomical and functional localizer to identify stimulus-driven voxels in V1 (Figure 2B). For those voxels that were significantly driven by the stimulus, we then calculated the phase lag, or relative response latency to the stimulus (Figure 2C). Consistent with prior studies (Birn et al., 2001; Lewis et al., 2016; de Zwart et al., 2005) we found a wide range of hemodynamic response lags within V1 (Figure 2A, bottom). From this distribution, we extracted groups of fast and slow responding voxels (Figure 2D–E). Then, for each voxel identified as fast or slow within the task run, we calculated that voxel’s frequency spectrum in the resting-state. Figure 2F shows the resting-state spectrum of a representative fast and slow voxel from a single subject where the differences in the spectra, particularly the distinct slope under 0.2 Hz, are visible.
Our simulations had predicted a difference in the overall frequency content of resting-state fMRI signals in voxels with fast versus slow HRFs. To quantify this property across voxels, we sought to generate a set of spectral features that could capture these resting-state spectral dynamics. We constructed four spectral features to capture spectral properties: the slope using a linear fit under 0.2 Hz, the exponent of an aperiodic 1 /f fit, the amplitude of low frequency fluctuations (0.01–0.1 Hz power; ALFF) (Zang et al., 2007), and the fractional ALFF (ratio of 0.01–0.08 Hz to 0–0.25 Hz; fALFF) (Zou et al., 2008; Figure 3). Each feature of the resting-state spectra revealed significant differences (Wilcoxon rank-sum test, p<0.05) between fast and slow voxels across subjects. The slope showed significant differences between the fast and slow voxels within each individual subject (15/15), while the aperiodic exponent, ALFF, and fALFF showed significant differences in 14, 11, and 13 subjects, respectively. (See Supplementary file 1 for p-values). Each of these features reflects information about the frequency content at high and low frequencies, suggesting that this was an effective metric for differentiating voxels with fast or slow hemodynamic responses. Notably, the most effective features for distinguishing fast and slow voxels were the ones that explicitly captured the relative difference in high-frequency vs. low-frequency power.
Faster hemodynamic responses in thalamus are also reflected in shallower frequency spectra
A key strength of fMRI is its ability to noninvasively image activity in the subcortex. Having established that the resting-state spectrum contained signatures of local hemodynamic response timing within V1, we next aimed to test whether this principle extended to the visual thalamus, specifically the LGN. Prior studies have shown that LGN has faster hemodynamic responses than V1 (Lau et al., 2011; Yen et al., 2011; Lewis et al., 2018); however, due to its small size and lower signal-to-noise ratio, extracting its spectral features accurately could be more challenging. We generated individual masks of the LGN within each subject using the individual anatomical segmentation (Fischl, 2012) and a functional localizer (Figure 4A). To confirm the presence of stimulus-locked oscillations in the LGN and to assess the relative timing of its response, we first examined the latency of the average response to the visual stimulus in the LGN, as compared to the fast and slow groups in V1. We observed that the LGN peaked before both the fast and slow groups in V1 (Figure 4B), consistent with prior work demonstrating faster hemodynamics in thalamus (Lau et al., 2011; Yen et al., 2011; Lewis et al., 2018). Then, to test whether this faster hemodynamic response was similarly linked to flatter resting-state spectra, we compared the LGN voxels’ resting-state features to the previously identified fast and slow voxels in the cortex. We found that for all subjects there were clear differences in each spectral feature within the LGN compared to both the fast and slow voxels of the visual cortex (Figure 4C–F; See Supplementary file 1 for p-values). The feature that had the poorest sensitivity to differences between the LGN features and cortical features is ALFF, which could be explained by ALFF’s higher sensitivity to non-neural noise sources (Zou et al., 2008). We thus observed even shallower frequency slopes for the fast LGN voxels – again consistent with our simulation results, demonstrating that this pattern held not just within V1 but even extended to the LGN of the thalamus. This observation was also robust to controlling for the higher thermal noise in LGN signals (Figure 4—figure supplement 1).
Resting-state spectral information better characterizes neurovascular coupling delays than a breath hold task
Perhaps the most established method of mapping hemodynamic latencies across the brain is using a breath hold task to quantify cerebrovascular reactivity (Chang et al., 2008; Kastrup et al., 1999; Liu et al., 2017; Bright et al., 2009; Chen et al., 2021; Pinto et al., 2021). Therefore, we tested whether similar information about hemodynamic latencies found using the resting-state spectra could be found in data from the breath hold task. To determine this, we first mapped the vascular latency in response to the breath hold task on a voxel-wise basis across the brain (Chang et al., 2008). We then compared each voxel’s relative vascular latency from the breath hold task across the three groups identified in the visual task – the fast and slow cortical voxels and LGN voxels. We found that some subjects showed the expected temporal sequence of activation, where LGN voxels and fast cortical voxels respond earlier than slow cortical voxels (Figure 5A), but this effect was not present in all subjects (Figure 5B). Even among the subjects that exhibited the expected order of activation, few of them showed individual-level significant differences in breath hold latency between the groups (Figure 5C). Specifically, in only 7/15 subjects was the average breath hold latency of fast cortical voxels significantly faster compared to slow cortical voxels. Furthermore, the breath hold latency in LGN voxels was slower than expected in some subjects: it was significantly slower than fast cortical voxels in 5 subjects and significantly slower than slow cortical voxels in 1 subject (Wilcoxon rank-sum test, p<0.05). These results demonstrate that although the vascular latencies derived from the breath hold task do show significant differences between fast and slow cortical voxels and LGN voxels across the group, this effect is less robust in individual subjects than the resting-state spectral features. Additionally, not all subjects demonstrated the expected order of latencies between the three groups – with LGN first followed by fast cortical and, lastly, slow cortical voxels. Taken together, these results suggest that the features of the resting-state spectrum capture additional information about local differences in neurovascular coupling delays.
Features of the resting-state spectrum can predict voxels with fast or slow hemodynamic response timing
Our results showed consistent signatures of hemodynamic response latency in the resting-state fMRI signal, suggesting that this information could potentially be used to predict local neurovascular latencies. Indeed, we found that these spectral features were significantly correlated with the absolute timing of their hemodynamic responses (Figure 3—figure supplement 3, Supplementary file 3). To investigate the utility of the resting-state spectrum to predict the temporal dynamics of the HRF, we tested whether support vector machines (SVMs) could classify slow, fast, and fastest (LGN) voxels using information from the resting-state spectrum. First, we trained a SVM to classify slow, fast, and LGN voxels based on the four features of the resting-state spectra identified in Figure 3. We found that the classifier validation accuracies, both within each individual subject and on the dataset that combined all subjects, were well above chance (Figure 6), demonstrating robust prediction of local hemodynamic delays (per-subject accuracy in Supplementary file 2). We next considered whether correlated information between neighboring voxels could be contributing this prediction, to test whether we could generalize to distant voxels. To control for instances of voxels in the test set being in close proximity to voxels in the training set, we trained a new set of models using voxels in a single hemisphere (left or right) as the training set, and voxels in the other hemisphere as the testing set. Following this procedure, the number of voxels for training and testing was often closer to a 50–50 split than an 80–20 split; however, validation accuracies both within individual subjects and across subjects nevertheless remained well above chance. Interestingly, the overall performance across subjects was higher, suggesting better generalization across subjects when training on spatially separated voxels (Figure 6—figure supplement 1).
Regression models aiming to predict specific response timing (rather than classifying fast vs. slow) were significantly above chance but still performed poorly, with an average coefficient of determination of 0.104 and an average root mean squared error of 0.923 (Supplementary file 4), likely due to the fact that absolute HRF timing depends on task state (Chen and Glover, 2015). We next investigated whether our chosen features were sufficient for classification, or whether additional useful information was present in the full resting-state spectrum, by training a second SVM classifier to use each voxel’s resting-state spectrum as the predictor for the classifier. The input to the classifier was the resting-state spectrum limited to up to 0.5 Hz to reduce the number of features fed into the model to avoid overfitting. Once again, we found that the accuracies, both within individual subjects and on the dataset that combined all subjects, were well above chance (Figure 6; per-subject accuracy in Supplementary file 2). The SVM classifiers trained using the spectrum trended toward slightly higher classification accuracies than the pre-selected features (Figure 6), but the performance of the two models was not significantly different (p=0.2828; Wilcoxon rank-sum test). We therefore concluded that our constructed spectral features are capturing key aspects of the resting-state spectrum that are important for predicting hemodynamic latency.
We next tested whether the resting-state signals performed better than the breath hold task at differentiating between fast, slow, and LGN voxels. We trained another SVM classifier to use the breath hold latencies as predictors and found that this model performed significantly above chance (Figure 6), consistent with prior work demonstrating its utility in predicting delays. However, it nevertheless performed significantly worse than both resting-state models (Wilcoxon rank-sum test, p<0.0005; per-subject accuracy in Supplementary file 2). This result further supported the conclusion that resting-state spectral information performs well at capturing differences in hemodynamic response timings.
Discussion
We conclude that the frequency spectrum of resting-state fMRI signals contains rich information about local hemodynamic timing. Our simulations first illustrated the rationale for this effect: HRFs with faster temporal dynamics produce less power in the low frequency bands and a shallower attenuation of power at high frequencies. This effect is preserved even when we account for the higher amplitude of slower HRFs and the 1 /f decay of the amplitude of spontaneous neural activity. Following this observation, we constructed several quantitative features based on the resting-state spectrum, each capturing a similar property: the relative amplitude of high vs. low frequencies, and when used as predictors for a SVM were able to classify individual voxels as fast or slow. Our findings demonstrate that the temporal properties of the HRF affect the spectral features of resting-state fMRI signals and present a framework for characterizing the temporal properties of the hemodynamic response across voxels, which is crucial for accurate fMRI analyses.
The work presented here has broad implications for fMRI studies using the frequency content of the BOLD signal to make inferences about intrinsic brain activity. Prior work has identified changes in the spectral content of fMRI signals in certain clinical populations, including major depressive disorder (Wang et al., 2016), mild cognitive impairment (Han et al., 2011), and Alzheimer’s disease (Yang et al., 2020), and interpreted these differences as changes in intrinsic brain activity. However, as seen here, differences in hemodynamic response timings will also alter the frequency content of fMRI signals. Furthermore, prior work has shown that intrinsic timescales of BOLD activity vary across brain regions and can be used to predict individual subject patterns (Fallon et al., 2020). Our work suggests hemodynamic differences may contribute to these observations. Therefore, changes in the spectral content of fMRI signals can arise not just from differences in intrinsic brain activity but could also be indicators of different hemodynamic response timing.
A common theme across the spectral features we constructed was that they were sensitive to the relative contributions of low- and high-frequency power. Two features we selected (the slope of the spectrum and the exponent of the aperiodic 1 /f fit) are direct measures of the attenuation in power towards higher frequencies. Similar information is contained in the commonly used metrics ALFF and fALFF (Zou et al., 2008): ALFF is a marker of the power in low frequency bands and fALFF is a measure of the ratio of high to low frequency power. While each of these features showed significant differences between fast and slow cortical voxels, only the slope and fALFF exhibited significant differences within each individual subject. These features explicitly capture the relative difference in high-frequency vs. low-frequency power, corresponding to the main prediction of our simulations. Conversely, ALFF had the poorest sensitivity within subjects which could be due to the fact that ALFF is the only feature that does not include information from both low and high frequency ranges and is instead limited to a relatively narrow band of frequencies, 0.01–0.08 Hz. Overall, our results suggest that while some information is contained in the magnitude of the fMRI signal alone, capturing the relative power in high versus low frequencies is the most important metric for predicting hemodynamic timing. Moreover, using the resting-state spectrum itself did not perform significantly better than our chosen features in classifying fast and slow voxels, suggesting that the relative power is indeed the primary contributor when predicting variations in local HRF timings.
Our results highlight the substantial variability in hemodynamic timing even within a single cortical area, and provide a way to predict this speed, but the specific reasons why individual voxels are fast or slow are not yet clear. A multitude of neurovascular and anatomical factors contribute to the variability in hemodynamic response timing across the brain. For example, the anatomical organization of the vascular network inherently introduces temporal dispersion to the hemodynamic response. As blood travels through the vascular network in response to changes in neuronal activity it not only encounters a range of vessel diameters, but also many branch points which each introduce delays. Prior work has shown that the timing and duration of the hemodynamic response both increase towards the cortical surface while the fastest, narrowest hemodynamic responses occurred in deeper gray matter (Siero et al., 2011; de Zwart et al., 2005; Duvernoy et al., 1981). This suggests the possibility that the fast responding voxels identified in our study are more concentrated in deeper gray matter, although the spatial resolution of our images was not sufficient to separate voxels by cortical depth in this study. Vascular elasticity and reactivity also change across the vascular network and vary between regions further contributing to the variable timing of the hemodynamic response across the brain (Uludağ and Blinder, 2018; Drew, 2019). Arteries are surrounded by smooth muscle cells that rapidly relax and expand their diameter in response to neural activity while veins dilate much more slowly (Drew, 2019). Local differences in the relative concentrations of arteries and veins could therefore also introduce differences in the temporal dynamics of the BOLD signal (Siero et al., 2011; de Zwart et al., 2005; Taylor et al., 2018), suggesting that the slow responding voxels may be more likely to contain veins. These relationships between vascular anatomy and timing further highlight that fMRI studies analyzing spectral power or temporal autocorrelation should take into consideration how local vascular anatomy affects these metrics.
A unique advantage of fMRI is its ability to image throughout the whole brain, and our results suggest the potential for extracting more information from fMRI studies of subcortex. Due to significant improvements in both the sensitivity and spatial resolution of fMRI, an increasing number of studies are utilizing fMRI to study small, deep brain structures such as the thalamus and brainstem (Beissner et al., 2014; Sclocco et al., 2018; Saranathan et al., 2021; Beissner and Baudrexel, 2014). The ability to image these deeper brain structures in humans opens the door to studying diverse aspects of cognition associated with these deeper brain regions (Beissner et al., 2014; Sclocco et al., 2018; Saranathan et al., 2021; Beissner and Baudrexel, 2014; Sherman, 2007). However, these deep brain structures also have unique physiological and anatomical properties that alter their vascular dynamics (Lau et al., 2011; Yen et al., 2011; Lewis et al., 2018; Duvernoy, 2009; Devonshire et al., 2012). Faster hemodynamic responses are frequently reported in subcortex, but hemodynamic responses in these regions are less well characterized than in the cortex. Our analyses found that variations in HRF timing are reflected in the frequency spectra of thalamic voxels as well, despite their differing signal-to-noise ratio. This result demonstrates the utility of our approach in subcortex and could benefit neuroimaging studies of structures such as the thalamus, a target of increasing interest in fMRI.
Improvements in the temporal resolution of BOLD fMRI have also sparked interest in detecting neural sequences at sub-second timescales, which are highly relevant for many studies of cognition. Recent studies have leveraged fast fMRI to detect rapid sequences of neural events related to visual sequence detection (Wittkuhn and Schuck, 2021), auditory dynamics (Frühholz et al., 2020), and changes in arousal state (Setzer et al., 2021). As we continue to identify these rapid neural sequences, it will become even more crucial to consider how the hemodynamic response varies across regions to determine whether a given sequence represents regional differences in neuronal or in hemodynamic timing. Considering spectral signatures can support inference of precise timing of neural activity by providing information about relative hemodynamic latencies between voxels and regions.
We found that the resting-state signals predicted voxel-wise differences in relative hemodynamic response timing significantly better using a breath hold to estimate vascular delays, the current gold standard. This observation could be explained by multiple factors. One potential factor is subject motion (Figure 5—figure supplement 1). There was significantly more motion during the breath hold task compared to the visual stimulus (Wilcoxon rank-sum test, p=0.01), which could lead to breath tasks providing less reliable information. Additionally, the models trained using the breath hold latencies were only able to leverage one piece of information: the measure of vascular latency derived from cross-correlation. The information from taking the cross-correlation in the time domain may also be more sensitive to noise, especially at the fast sampling rate of our scans. However, a more important factor may be the biological mechanisms generating each signal. The breath hold task directly modulates cerebrovascular reactivity with minimal accompanying changes in CMRO2, allowing for assessment of local cerebral vascular reactivity uncoupled from neuronal activation (Kastrup et al., 1999). CVR is important to assess in many clinical applications, and the breath hold-based approach remains a gold standard for CVR mapping (Kastrup et al., 1999; Liu et al., 2017; Bright et al., 2009; Pinto et al., 2021). However, while CVR is a significant modulator of the hemodynamic response, it is only one component of neurovascular coupling, and there are many other factors and signaling pathways that also contribute. In particular, local metabolic factors and feedback mechanisms also modulate of blood flow in the brain and are not replicated in the breath hold task (Iadecola, 2017; Hosford and Gourine, 2019). The hemodynamic responses induced by neural activity may, therefore, not be identical in timing to those induced by a purely vascular signal. By contrast, the signals obtained from resting-state fMRI are coupled to underlying neural activity in an analogous manner as during a task condition (Mateo et al., 2017; Ma et al., 2016; Fox and Raichle, 2007). The improved performance of the resting-state-based prediction may therefore reflect that it is intrinsically more similar to a task than the breath hold condition, as the underlying biological origins of resting-state signals share common mechanisms with task-induced neurovascular coupling. Future applications may therefore benefit from continuing to use breath tasks as a gold standard to assess CVR, whereas resting-state analyses may be a better metric of neurovascular coupling.
One limitation of our study is that our approach has thus far only been validated in the visual system, which is the basis for the majority of our knowledge about the HRF. Given that the HRF has been shown to vary across the brain, the generalizability of our approach to different brain regions has yet to be established. However, studies that have empirically derived the HRF in other brain regions, such as the motor cortex (Siero et al., 2011) or somatosensory cortex (de Zwart et al., 2005), have produced HRFs with similar shapes and parameters to the double-gamma HRF utilized in our simulations (Handwerker et al., 2004; Siero et al., 2011; de Zwart et al., 2005; Taylor et al., 2018) providing evidence that this approach may be generalizable across the brain. Still, future work could investigate this approach in other primary sensory systems such as the auditory system where there are reliable stimulation paradigms. Additionally, this study used data collected at an ultra-high magnetic field strength (7T) which affords a higher signal-to-noise ratio (SNR) compared to more traditional 3T fMRI. However, we also performed the same spectral analysis in an independent dataset acquired at 3T and robustly replicated our result, demonstrating that this pattern is preserved despite a lower SNR and different acquisition parameters (Figure 3—figure supplement 1).
The work discussed here suggests a wide range of neuroscience applications for our approach to measuring hemodynamic timing. One logical next step would be to use the ability to characterize temporal variation in the HRF to not just predict, but to correct for vascular delays. Previous work demonstrated that correcting for varying hemodynamic latencies across the brain can affect functional connectivity analyses (Rangaprakash et al., 2018; Chang et al., 2008). Our results could further enhance removal of non-neural latency differences that confound functional connectivity metrics, both static and dynamic, and can increase confidence that the networks we are analyzing are derived from neuronal dynamics. This has become increasingly of interest as more studies use functional connectivity, particularly in resting-state, to study dynamics underlying diseases such as PTSD (Jin et al., 2017; Rabinak et al., 2011; Maron-Katz et al., 2020; Zhu et al., 2017), Alzheimer’s Disease (Greicius et al., 2004; Bai et al., 2008; Agosta et al., 2012; Koch et al., 2012), Parkinson’s Disease (Helmich et al., 2010; Hacker et al., 2012; Agosta et al., 2014; Baggio et al., 2015; Putcha et al., 2015), and others (Filippi et al., 2019). Since changes in neurovascular function have been observed in many disorders (de la Torre and Mecocci, 2018; Gutteridge et al., 2020; Burrage et al., 2018), analyzing spectral dynamics may help interpret functional connectivity differences in clinical populations.
Although we focused on the utility of resting-state spectral information to classify voxels as having relatively fast or slow hemodynamic response timings, we did also examine the correlation with absolute timing measures. We found that most subjects had significant positive correlations between each of the spectral features and the task response latency. Despite this fact, when estimating continuous predictions of absolute hemodynamic response lag, prediction performance was poor. Importantly, previous evidence has shown that the absolute timing of the HRF varies between task and resting-state conditions (Chen and Glover, 2015), suggesting that caution is needed in using absolute measures of hemodynamic latencies, as these may not generalize across conditions. Measuring relative differences in HRF timing might therefore be a more effective method to correct for variations in the hemodynamic response across brain regions, since task state can modulate the absolute timing of the hemodynamic response. Furthermore, the relative differences provide the key information necessary to interpret sequences of activity across brain regions, enabling examination of whether purely vascular differences are present across those regions.
Together, our results demonstrate that the resting-state fMRI signal contains information about local hemodynamic response speeds. This approach can help understand brain-wide variations in HRF dynamics, which is critical as the field moves toward a new era of fMRI studies utilizing fast fMRI to study rapid neuronal dynamics and higher level cognition.
Materials and methods
Simulations
Request a detailed protocolSpectra of simulated BOLD responses were generated by convolving a given HRF with oscillating stimuli, ranging from 0.1 to 0.5 Hz, and taking the magnitude of the simulated BOLD response as the power at that frequency (Figure 1A). We used six HRFs with varying time-to-peak (TTP), full width at half maximum (FWHM), and peak percent signal changes (PSCs) to represent a range of physiologically relevant HRFs (Figure 1B). These properties were drawn from previous work characterizing varying HRF temporal dynamics at different cortical depths (Siero et al., 2011) and the values used are reported below in Supplementary file 5. We also normalized these HRFs by their maximum percent signal change and re-simulated the BOLD responses to create a new simulated spectrum for each HRF (Figure 1C). Additionally, we performed simulations to account for the dominant 1/f-like spectral pattern of neural activity by setting the amplitude of the oscillating stimuli to be 1/stimulus frequency ranging from 0.1 to 0.5 Hz (Figure 1—figure supplement 1).
Subject population
Request a detailed protocolAll experimental procedures were approved by the Massachusetts General Hospital Institutional Review Board (protocol number: 2014P001068), and all subjects provided informed consent. Twenty-one participants were scanned in total; five were excluded for excessive motion and one was excluded for poor performance on the visual task suggesting they had closed their eyes. This left 15 subjects whose data was analyzed (mean age = 28 years, range = 22–42 years, 8 female).
Experimental design
Request a detailed protocolSubjects underwent a total of 7 functional scans: 3 visual stimulus, 2 breath hold, and 2 resting-state runs. All stimuli were programmed in MATLAB using Psychtoolbox (Brainard, 1997).
Visual stimulus
Request a detailed protocolEach visual stimulus functional run lasted 254 s, with the first 14 s showing a gray screen with the red fixation dot and the following 240 s consisting of the 12 Hz counterphase flickering radial checkerboard. To drive continuous neural oscillations in the visual cortex, the luminance contrast of the flickering checkerboard oscillated at a frequency of 0.05 Hz (except for one subject who was presented with 0.1 Hz oscillations). To assist the subjects with fixation, in the center of the visual field was a red dot that changed brightness at random intervals. Subjects were directed to press a button whenever the brightness of the red dot changed, and their average response time and response accuracy was reported at the end of each run. This allowed us to monitor participant engagement with the task. Each subject participated in 3 visual stimulus runs.
Breath hold Task
Request a detailed protocolFor each breath hold run subjects performed eight repetitions of an adapted version of a previously established breath hold task (Chang et al., 2008): a block comprised of 27 s of free breathing, 3 cycles of paced breathing (3 s breathe out, 3 s breathe in), a 15 s breath hold, and, lastly, a 30-s period of free breathing. The total time for the breath hold task scans was 8.5 min. The instructions for breathing were projected to the subjects displaying the text ‘Breathe Freely’; ‘Breathe Out’; ‘Breathe In’; and ‘HOLD BREATH’. All but one subject participated in 2 breath hold runs; a single subject performed only 1 breath hold run. For 2 subjects, a single breath hold run was excluded from analyses for excessive motion defined here as greater than 0.5 mm average motion across the whole run. No additional runs were excluded for these subjects.
Resting-state
Request a detailed protocolEach subject participated in 2 resting-state runs where they were instructed to relax in the scanner with their eyes open and to try not to fall asleep. For some of the subjects, a fixation dot was presented to help minimize eye movements. Each resting-state scan lasted 8.5 minutes.
MRI Data Acquisition
Request a detailed protocolSubjects were scanned on a 7 Tesla Siemens MAGNETOM scanner with a custom-built 32-channel head coil. Anatomical images were acquired with 0.75 mm isotropic multiecho magnetization-prepared rapid gradient-echo (MEMPRAGE) protocol (van der Kouwe et al., 2008) with TR = 2.530ms, echo time (TE)=1.76ms and 3.7ms, inversion time (TI)=1100ms, echo-spacing=6.2ms, 7o flip angle, bandwidth = 651 Hz, in-plane acceleration R=2, FOV = 320 x 320 x 244 mm and a total scan time of 7:20 min. For functional runs, 15 oblique slices were positioned to target the calcarine sulcus to include primary visual cortex (V1) and angled to include the lateral geniculate nucleus (LGN) located in the thalamus. Functional runs were acquired as single-shot gradient-echo EPI with 2 mm isotropic resolution, TR = 227ms, TE = 24ms, echo-spacing=0.59ms, 30o flip angle, bandwidth = 2604 Hz, in-plane acceleration R=2, SMS Multiband Factor = 3, CAIPI shift = FOV/3 (Setsompop et al., 2012).
Physiological Monitoring
Request a detailed protocolFor all the functional scans, subjects’ heart rate and respiration were monitored using piezoelectric transducer on the non-dominant thumb and a respiratory belt around the upper rib cage, respectively. The physiological recordings were obtained at a sampling rate of 1000 Hz using a PowerLab physio box connected to a computer running LabChart 7 from ADInstruments.
fMRI Analyses
fMRI preprocessing
Request a detailed protocolAnatomical images were bias-corrected using SPM (https://www.fil.ion.ucl.ac.uk/spm/) and segmented using FreeSurfer (Fischl, 2012). Functional runs were preprocessed with slice-timing correction, performed using FSL (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/), and motion correction, performed using AFNI software (https://afni.nimh.nih.gov/). No spatial smoothing was applied.
Because fast fMRI has distinct contributions from systemic physiological noise, including cardiac rhythms and respiration, physiological noise removal was performed on the visual stimulus and resting-state functional runs using dynamic regression adapted from RETROICOR (Glover et al., 2000) in runs where physiological recordings were successfully collected. In runs where physiological recordings were not successfully collected (2 runs total), physiological noise was removed using a statistical model of harmonic regression with autoregressive noise (HRAN) (Agrawal et al., 2020).
Visual localizer
Request a detailed protocolFor each subject, one of the visual stimulus runs was used as a functional localizer to identify voxels that were significantly driven by the oscillating stimulus. A general linear model (GLM) was fit in FSL (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/) using sine and cosine basis functions with the same period as the stimulus. The F-statistic of the combined fit to both the sine and cosine basis function was transformed to a Z-score and voxels with a Z-score above 2.5 were selected for further analysis. This functional localizer was then constrained by its intersection with the anatomical definition of V1 (Figure 2B). Specifically, the V1 segmentation was generated automatically from the MEMPRAGE volume based on the cortical surface reconstruction generated using FreeSurfer (Fischl, 2012). The selected voxels from the localizer run were then mapped to each other functional run in a single transformation step. This was done by first registering all functional runs to the anatomical scan using boundary-based registration (Greve and Fischl, 2009) and then resampling the desired volume into the localizer field of view using the registration matrices.
LGN segmentation
Request a detailed protocolThe lateral geniculate nucleus (LGN) was segmented using both anatomical and functional constraints. The anatomically defined boundaries of the LGN were generated using the FreeSurfer developmental version that generates an individual-level probabilistic atlas in individual anatomical space (Iglesias et al., 2018). From this probabilistic atlas we considered voxels with at least a 30% probability of falling within LGN and dilated this mask to capture border voxels. We next applied a functional constraint using the visual localizer where voxels in the dilated mask with a Z-score above 2.5 were considered part of our final LGN map in each subject.
Voxel-wise phase analysis and groupings
Request a detailed protocolWe then averaged the two remaining visual stimulus runs and extracted voxel’s average time series between the two runs. We discarded the first 14 s to analyze the steady-state response to the visual stimulus. An estimation of the voxel’s lag in relation to the oscillating stimulus was calculated using the arctangent of the sine and cosine regressor estimates. This allowed us to generate a histogram of latencies to the visual stimulus (Figure 2C). To extract ‘fast’ and ‘slow’ reacting voxels, a Gaussian model was fit to the histogram of phase delays. The centroid () and Full Width at Half Maximum () of the Gaussian fit were calculated. Groups of fast and slow voxels were then Edges of the Gaussian fit were defined as ± and fast and slow groups were made that each had a width of . Fast voxels were identified as being within while slow voxels were identified as (Figure 2C). This procedure was done for each individual subject and on average yielded 142 fast voxels (range 64–236) and 139 slow voxels (range 68–254). Masks of fast and slow voxels were generated per subject and then mapped to the resting-state runs to inform the spectral analysis (Figure 2D).
Resting-state spectral analysis
Fast and slow voxels were always identified in task-driven runs, allowing us to assess frequency content in the resting-state run using fully independent data. The maps of fast and slow voxels were registered to each individual resting-state run, and for each voxel within these masks, after discarding the first 14 s, the voxel-wise resting-state power spectrum was calculated using the Chronux toolbox (Bokil et al., 2010) with five tapers. We used four features to characterize the resting-state spectra: (1) slope of linear fit under 0.2 Hz; (2) the exponent of the aperiodic 1 /f fit under 0.5 Hz; (3) the amplitude of low frequency fluctuations (ALFF); and (4) the fractional ALFF (fALFF). Each of these features was z-scored within the run and then averaged on a voxel-wise basis between the two resting-state runs. All analysis of the resting-state spectra was performed in MATLAB. See Figure 3A–D for more information.
Slope
Request a detailed protocolA linear fit was generated for each voxel’s resting-state spectrum under 0.2 Hz using least-squares to determine the coefficients of a first order polynomial. From this we were able to record the slope of that linear fit for each voxel.
Exponent of aperiodic fit
Request a detailed protocolEquation 1 was fit using the Levenberg-Marquardt algorithm to solve non-linear least squares for each voxel’s resting-state spectra. In this equation, is the independent variable and the and are the values being fit. The exponent of the resultant fit () was recorded for each voxel.
Amplitude of low-frequency fluctuations (ALFF)
Request a detailed protocolFor each voxel, ALFF was calculated according to the method outlined in Zou et al., 2008. Briefly, each voxel’s time series was band pass filtered between 0.01 and 0.08 Hz. Then, the voxel’s time series is transformed into the frequency domain via Fast Fourier Transform (FFT). The power at each frequency is proportional to the square of the amplitude of the FFT at that frequency, and for each voxel, the ALFF value was taken as the averaged square root of the power in the 0.01–0.08 Hz frequency range. This is shown in Equation 2 where is the magnitude of the FFT at frequency , is the number of time points, is the power spectrum, and the bar represents the average in the specified range.
Fractional ALFF (fALFF)
Request a detailed protocolEach voxel’s fALFF was calculated as described in Zou et al., 2008. Fractional ALFF is briefly defined as the ratio of the power of each frequency at the low frequency range (0.01–0.08 Hz) to that of the ‘global’ frequency range (0.01–0.25 Hz). See Equation 3.
Breath-hold latency calculations
Request a detailed protocolVoxel-wise hemodynamic latencies were calculated according to the method outlined in Chang et al., 2008. Each voxel’s hemodynamic latency was defined as the time-lag yielding the maximum cross-correlation between the given voxel’s time series, , and a reference time series, . The reference time series was found by taking the average time series across voxels in the brain that exceeded a minimum correlation of with the breath hold task regressor. This breath hold task regressor was defined as the convolution of a box car function, where the value is set to 1 during the breath hold and 0 at other times, and a sign-reversed canonical HRF (Glover, 1999). Both and were resampled to a resolution of 100ms before computing cross-correlations.
Support vector machine (SVM) classification of fast, slow, and LGN voxels
Request a detailed protocolSVM classifiers were trained both within and across subjects using Scikit-learn in Python (Pedregosa, 2011). Three models were trained with different predictive features: (1) the resting-state spectrum between 0 and 0.5 Hz, (2) the 4 features of the resting-state spectrum previously identified, and (3) the latency of the response to the breath hold task. The resting-state spectrum was extracted by taking the power at frequencies up to 0.5 Hz ultimately generating a set of 461 features. For all models, before being put into the classifier, the data was normalized by removing the mean and dividing by the standard deviation across voxels for each feature independently using the StandardScaler function of Scikit-learn. For all models, the parameters of the SVM classifier were as follows: regularization parameter ()=10, kernel type =radial basis function (rbf), kernel coefficient . To get validation accuracies both within and across subjects, 1000 bootstraps were performed where the 80–20 test-train split of voxels was randomly chosen for each bootstrap. The average validation accuracies over the 1000 bootstraps were calculated along with the 95% confidence intervals. This methodology was followed for all 3 SVM classifiers.
Support vector machine (SVM) regression to predict relative hemodynamic response latency
Request a detailed protocolSVM models for regression were trained within subject to continuously predict the hemodynamic response latency estimated from the visual stimulus. We limited the regression to voxels whose hemodynamic response latency relative to the median was between [–3, 3] sec. Two different models were trained with different predictive features: (1) the resting-state spectrum between 0 and 0.5 Hz and (2) the 4 features of the resting-state spectrum previously identified. The input data to the model was normalized by removing the mean and dividing by the standard deviation cross voxels using the StandardScaler function of Scikit-learn. The parameters of the SVM regression model were as follows: kernel type = radial basis function, kernel coefficient , regularization parameter = 10, epsilon = 0.1. To get the validation coefficient of determination (, 1000 bootstraps were performed on an 80–20 test-train split of voxels that were randomly chosen for each bootstrap. The average over the 1000 bootstraps was calculated along with the 95% confidence interval.
Replication dataset at 3 Tesla
Request a detailed protocolTo generate the replication figure (Figure 3—figure supplement 1), we analyzed 10 subjects (mean age = 23 years, range = 19–29, 6 female) from a previously published dataset (Williams et al., 2023). Subjects were selected for analysis if they had completed at least one run of a visual task as well as one resting-state scan of at least 10 min. Briefly, subjects were scanned on a 3 Tesla Siemens Prisma scanner with a 64-channel head and neck coil. Anatomical runs were acquired with 1 mm isotropic T1-weighted multi-echo MPRAGE (van der Kouwe et al., 2008). Functional runs were acquired using TR = 0.378 s, 2.5 mm isotropic voxels, and multiband factor = 8. Similar preprocessing steps were performed for this dataset as described above in fMRI Preprocessing including slice-timing correction and motion correction. More information about the specific acquisition parameters and preprocessing can be found in the original open-access publication (Williams et al., 2023).
The visual stimulus was a counterphase 12 Hz flickering checkerboard stimulus that lasted 254 s with fixed 16 s ON and 16 s OFF periods beginning with an OFF period. To assist the subjects with fixation, in the center of the checkerboard was a red dot that changed brightness at random interval. Subjects were instructed to press a button whenever they detected a color change. During the resting-state scans subjects performed a behavioral task where they were instructed to close their eyes and press a button on a response box on every breath in. This behavioral task enabled us to behaviorally monitor if the subjects fell asleep. Only subjects that were consistently doing the behavioral task during the entire resting-state run were included in our analysis.
A GLM was fit in FSL to identify voxels that were driven by the stimulus and voxels with a Z-statistic above 2.5 were selected for further analysis. This functional localizer was then further constrained by its intersection with the anatomical definition of V1. Voxel-wise hemodynamic response lags within this combined mask were calculated and groups of fast and slow voxels were extracted using the methods described above in Voxel-wise Phase Analysis and Groupings. Then, the fast and slow voxels were registered to the resting-state run and the same procedure outlined above in Resting-state Spectral Analysis was used to extract the 4 spectral features of interest.
Data availability
The data used in this paper has been deposited on OpenNeuro (https://doi.org/10.18112/openneuro.ds004645.v1.0.0).
-
OpenNeuroResting-state fMRI signals contain spectral signatures of local hemodynamic response timing.https://doi.org/10.18112/openneuro.ds004645.v1.0.0
References
-
Resting state fMRI in Alzheimer’s disease: beyond the default mode networkNeurobiology of Aging 33:1564–1578.https://doi.org/10.1016/j.neurobiolaging.2011.06.007
-
Cortico-striatal-thalamic network functional connectivity in hemiparkinsonismNeurobiology of Aging 35:2592–2602.https://doi.org/10.1016/j.neurobiolaging.2014.05.032
-
Cognitive impairment and resting-state network connectivity in Parkinson’s diseaseHuman Brain Mapping 36:199–212.https://doi.org/10.1002/hbm.22622
-
Simultaneous multislice (SMS) imaging techniquesMagnetic Resonance in Medicine 75:63–81.https://doi.org/10.1002/mrm.25897
-
Investigating the human brainstem with structural and functional MRIFrontiers in Human Neuroscience 8:116.https://doi.org/10.3389/fnhum.2014.00116
-
Chronux: A platform for analyzing neural signalsJournal of Neuroscience Methods 192:146–151.https://doi.org/10.1016/j.jneumeth.2010.06.020
-
Cerebrovascular dysfunction with stress and depressionBrain Circulation 4:43–53.https://doi.org/10.4103/bc.bc_6_18
-
Altered cerebrovascular reactivity due to respiratory rate and breath holding: a BOLD-fMRI study on healthy adultsBrain Structure & Function 226:1229–1239.https://doi.org/10.1007/s00429-021-02236-5
-
The vascular hypothesis of alzheimer’s disease: a key to preclinical prediction of dementia using neuroimagingJournal of Alzheimer’s Disease 63:35–52.https://doi.org/10.3233/JAD-180004
-
Alterations in the BOLD fMRI signal with ageing and disease: a challenge for neuroimagingNature Reviews. Neuroscience 4:863–872.https://doi.org/10.1038/nrn1246
-
Vascular and neural basis of the BOLD signalCurrent Opinion in Neurobiology 58:61–69.https://doi.org/10.1016/j.conb.2019.06.004
-
Cortical blood vessels of the human brainBrain Research Bulletin 7:519–579.https://doi.org/10.1016/0361-9230(81)90007-1
-
Spontaneous fluctuations in brain activity observed with functional magnetic resonance imagingNature Reviews. Neuroscience 8:700–711.https://doi.org/10.1038/nrn2201
-
BOLD neurovascular coupling does not change significantly with normal agingHuman Brain Mapping 38:3538–3551.https://doi.org/10.1002/hbm.23608
-
Cerebrovascular function during cognition in Parkinson’s disease: A functional transcranial Doppler sonography studyJournal of the Neurological Sciences 408:116578.https://doi.org/10.1016/j.jns.2019.116578
-
Spatial remapping of cortico-striatal connectivity in Parkinson’s diseaseCerebral Cortex 20:1175–1186.https://doi.org/10.1093/cercor/bhp178
-
What is the key mediator of the neurovascular coupling response?Neuroscience and Biobehavioral Reviews 96:174–181.https://doi.org/10.1016/j.neubiorev.2018.11.011
-
Dynamic brain connectivity is a better predictor of PTSD than static connectivityHuman Brain Mapping 38:4479–4496.https://doi.org/10.1002/hbm.23676
-
Assessment of cerebral oxidative metabolism with breath holding and fMRIMagnetic Resonance in Medicine 42:608–611.https://doi.org/10.1002/(sici)1522-2594(199909)42:3<608::aid-mrm26>3.0.co;2-i
-
Cerebrovascular pathophysiology following mild traumatic brain injuryClinical Physiology and Functional Imaging 31:85–93.https://doi.org/10.1111/j.1475-097X.2010.00990.x
-
Individual Patterns of Abnormality in Resting-State Functional Connectivity Reveal Two Data-Driven PTSD SubgroupsAmerican Journal of Psychiatry 177:244–253.https://doi.org/10.1176/appi.ajp.2019.19010060
-
Scikit-learn: machine learning in pythonJournal of Machine Learning Research: JMLR 12:2825–2830.
-
Cerebrovascular reactivity mapping without gas challenges: a methodological guideFrontiers in Physiology 11:1711.https://doi.org/10.3389/fphys.2020.608475
-
Hemodynamic response function (HRF) variability confounds resting-state fMRI functional connectivityMagnetic Resonance in Medicine 80:1697–1713.https://doi.org/10.1002/mrm.27146
-
Blipped-controlled aliasing in parallel imaging for simultaneous multislice echo planar imaging with reduced g-factor penaltyMagnetic Resonance in Medicine 67:1210–1224.https://doi.org/10.1002/mrm.23097
-
Rapid brain MRI acquisition techniques at ultra-high fieldsNMR in Biomedicine 29:1198–1221.https://doi.org/10.1002/nbm.3478
-
The thalamus is more than just a relayCurrent Opinion in Neurobiology 17:417–422.https://doi.org/10.1016/j.conb.2007.07.003
-
Cortical depth-dependent temporal dynamics of the BOLD response in the human brainJournal of Cerebral Blood Flow and Metabolism 31:1999–2008.https://doi.org/10.1038/jcbfm.2011.57
-
Cortical depth dependence of the BOLD initial dip and poststimulus undershoot in human visual cortex at 7 TeslaMagnetic Resonance in Medicine 73:2283–2295.https://doi.org/10.1002/mrm.25349
-
Nonparametric hemodynamic deconvolution of FMRI using homomorphic filteringIEEE Transactions on Medical Imaging 34:1155–1163.https://doi.org/10.1109/TMI.2014.2379914
-
An improved approach to detection of amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI: fractional ALFFJournal of Neuroscience Methods 172:137–141.https://doi.org/10.1016/j.jneumeth.2008.04.012
Article and author information
Author details
Funding
National Institutes of Health (R00-MH111748)
- Laura D Lewis
National Institutes of Health (U19-NS123717)
- Laura D Lewis
National Institutes of Health (R01-AG070135)
- Laura D Lewis
Searle Scholars Program
- Laura D Lewis
Pew Charitable Trusts (Pew Biomedical Scholars Program)
- Laura D Lewis
Sloan Research Fellowship (Sloan Fellowship)
- Laura D Lewis
One Mind (One Mind Rising Star Award)
- Laura D Lewis
National Institutes of Health (T32-GM008764)
- Sydney M Bailes
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This work was supported by National Institutes of Health grants R00-MH111748, U19-NS123717, and R01-AG070135, the Searle Scholars Program, the Pew Biomedical Scholars Program, the Sloan Fellowship, and the One Mind Rising Star award. Resources were provided by NIH grant P41-EB030006 and T32-GM008764.
Ethics
All experimental procedures were approved by the Massachusetts General Hospital Institutional Review Board and all subjects provided informed consent. (protocol number: 2014P001068).
Copyright
© 2023, Bailes et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 1,263
- views
-
- 157
- downloads
-
- 6
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Neuroscience
Granule cells of the cerebellum make up to 175,000 excitatory synapses on a single Purkinje cell, encoding the wide variety of information from the mossy fibre inputs into the cerebellar cortex. The granule cell axon is made of an ascending portion and a long parallel fibre extending at right angles, an architecture suggesting that synapses formed by the two segments of the axon could encode different information. There are controversial indications that ascending axon (AA) and parallel fibre (PF) synapse properties and modalities of plasticity are different. We tested the hypothesis that AA and PF synapses encode different information, and that the association of these distinct inputs to Purkinje cells might be relevant to the circuit and trigger plasticity, similar to the coincident activation of PF and climbing fibre inputs. Here, by recording synaptic currents in Purkinje cells from either proximal or distal granule cells (mostly AA and PF synapses, respectively), we describe a new form of associative plasticity between these two distinct granule cell inputs. We show for the first time that synchronous AA and PF repetitive train stimulation, with inhibition intact, triggers long-term potentiation (LTP) at AA synapses specifically. Furthermore, the timing of the presentation of the two inputs controls the outcome of plasticity and induction requires NMDAR and mGluR1 activation. The long length of the PFs allows us to preferentially activate the two inputs independently, and despite a lack of morphological reconstruction of the connections, these observations reinforce the suggestion that AA and PF synapses have different coding capabilities and plasticity that is associative, enabling effective association of information transmitted via granule cells.
-
- Neuroscience
Rats are believed to communicate their emotional state by emitting two distinct types of ultrasonic vocalizations. The first is long ‘22-kHz’ vocalizations (>300 ms, <32-kHz) with constant frequency, signaling aversive states, and the second is short ‘50-kHz’ calls (<150 ms, >32 kHz), often frequency-modulated, in appetitive situations. Here, we describe aversive vocalizations emitted at a higher pitch by male Wistar and spontaneously hypertensive rats (SHR) in an intensified aversive state – prolonged fear conditioning. These calls, which we named ‘44-kHz’ vocalizations, are long (>150 ms), generally at a constant frequency (usually within 35–50-kHz range) and have an overall spectrographic image similar to 22-kHz calls. Some 44-kHz vocalizations are comprised of both 22-kHz-like and 44-kHz-like elements. Furthermore, two separate clustering methods confirmed that these 44-kHz calls can be separated from other vocalizations. We observed 44-kHz calls to be associated with freezing behavior during fear conditioning training, during which they constituted up to 19.4% of all calls and most of them appeared next to each other forming uniform groups of vocalizations (bouts). We also show that some of rats’ responses to the playback of 44-kHz calls were more akin to that of aversive calls, for example, heart rate changes, whereas other responses were at an intermediate level between aversive and appetitive calls. Our results suggest that rats have a wider vocal repertoire than previously believed, and current definitions of major call types may require reevaluation. We hope that future investigations of 44-kHz calls in rat models of human diseases will contribute to expanding our understanding and therapeutic strategies related to human psychiatric conditions.