Europe PMC

This website requires cookies, and the limited processing of your personal data in order to function. By using the site you are agreeing to this as outlined in our privacy notice and cookie policy.

Abstract 


The blood flow response to a vasoactive stimulus demonstrates regional heterogeneity across both the healthy brain and in cerebrovascular pathology. The timing of a regional hemodynamic response is emerging as an important biomarker of cerebrovascular dysfunction, as well as a confound within fMRI analyses. Previous research demonstrated that hemodynamic timing is more robustly characterized when a larger systemic vascular response is evoked by a breathing challenge, compared to when only spontaneous fluctuations in vascular physiology are present (i.e., in resting-state data). However, it is not clear whether hemodynamic delays in these two conditions are physiologically interchangeable, and how methodological signal-to-noise factors may limit their agreement. To address this, we generated whole-brain maps of hemodynamic delays in nine healthy adults. We assessed the agreement of voxel-wise gray matter (GM) hemodynamic delays between two conditions: resting-state and breath-holding. We found that delay values demonstrated poor agreement when considering all GM voxels, but increasingly greater agreement when limiting analyses to voxels showing strong correlation with the GM mean time-series. Voxels showing the strongest agreement with the GM mean time-series were primarily located near large venous vessels, however these voxels explain some, but not all, of the observed agreement in timing. Increasing the degree of spatial smoothing of the fMRI data enhanced the correlation between individual voxel time-series and the GM mean time-series. These results suggest that signal-to-noise factors may be limiting the accuracy of voxel-wise timing estimates and hence their agreement between the two data segments. In conclusion, caution must be taken when using voxel-wise delay estimates from resting-state and breathing-task data interchangeably, and additional work is needed to evaluate their relative sensitivity and specificity to aspects of vascular physiology and pathology.

Free full text 


Logo of nihpaLink to Publisher's site
Neuroimage. Author manuscript; available in PMC 2023 Jul 1.
Published in final edited form as:
PMCID: PMC10208394
NIHMSID: NIHMS1900796
PMID: 37072074

Hemodynamic timing in resting-state and breathing-task BOLD fMRI

Associated Data

Supplementary Materials
Data Availability Statement

Abstract

The blood flow response to a vasoactive stimulus demonstrates regional heterogeneity across both the healthy brain and in cerebrovascular pathology. The timing of a regional hemodynamic response is emerging as an important biomarker of cerebrovascular dysfunction, as well as a confound within fMRI analyses. Previous research demonstrated that hemodynamic timing is more robustly characterized when a larger systemic vascular response is evoked by a breathing challenge, compared to when only spontaneous fluctuations in vascular physiology are present (i.e., in resting-state data). However, it is not clear whether hemodynamic delays in these two conditions are physiologically interchangeable, and how methodological signal-to-noise factors may limit their agreement. To address this, we generated whole-brain maps of hemodynamic delays in nine healthy adults. We assessed the agreement of voxel-wise gray matter (GM) hemodynamic delays between two conditions: resting-state and breath-holding. We found that delay values demonstrated poor agreement when considering all GM voxels, but increasingly greater agreement when limiting analyses to voxels showing strong correlation with the GM mean time-series. Voxels showing the strongest agreement with the GM mean time-series were primarily located near large venous vessels, however these voxels explain some, but not all, of the observed agreement in timing. Increasing the degree of spatial smoothing of the fMRI data enhanced the correlation between individual voxel time-series and the GM mean time-series. These results suggest that signal-to-noise factors may be limiting the accuracy of voxel-wise timing estimates and hence their agreement between the two data segments. In conclusion, caution must be taken when using voxel-wise delay estimates from resting-state and breathing-task data interchangeably, and additional work is needed to evaluate their relative sensitivity and specificity to aspects of vascular physiology and pathology.

Keywords: BOLD fMRI, Hemodynamics, Resting-state, Breath-hold, Relative timing

1. Introduction

Oxygenated blood takes time to arrive at the brain, traveling from the heart and through the vascular tree to ultimately perfuse into local brain tissue where nutrient exchange can occur. Changes in blood flow can lead to changes in blood oxygenation, particularly in the venous vessels downstream of the exchange with tissue. Functional Magnetic Resonance Imaging (fMRI) is sensitive to blood flow changes due to their impact on local blood oxygenation levels, referred to as Blood Oxygenation Level Dependent (BOLD) contrast (Glover, 2011). However, processes which initiate local blood flow changes will not produce instantaneous changes in BOLD fMRI signals. A single neural event evokes a local hemodynamic response through neurovascular coupling mechanisms (Huneau et al., 2015), and this response is typically modeled with the canonical hemodynamic response function or HRF (Friston et al., 1998; Martindale et al., 2003; Penny et al., 2006). Determined empirically, this HRF peaks many seconds after an impulse of neural stimulation suggesting that there is an inherent BOLD fMRI “delay” due to the evolution of the local blood flow response to nearby neural activity. There is evidence that the shape and timing of the HRF to a neural event may vary across brain regions (Handwerker et al., 2004), with healthy aging (Issard and Gervain, 2018; West et al., 2019), and under different baseline states (Cohen et al., 2002; Liu et al., 2004). Via a different process to neurovascular coupling, blood flow changes can be evoked systemically by a non-neural vasodilatory stimulus such as carbon dioxide, a phenomenon known as cerebrovascular reactivity (CVR) (Chen and Gauthier, 2021; Liu et al., 2019; Pinto et al., 2020; Sleight et al., 2021). These CVR responses are also delayed, with regional heterogeneity in the vascular response observed across the healthy brain and altered by disease (Bright et al., 2009; Donahue et al., 2016; Holmes et al., 2020; Leung et al., 2016; Moia et al., 2020; Siegel et al., 2016; Sousa et al., 2014). There are several additional non-neural phenomena thought to drive variability in the spontaneous and systemic low frequency oscillation (sLFO) in BOLD contrast. Typically defined below 0.15 Hz, evidence of sources of variation in this non-stationary sLFO point towards blood CO2, heart rate and respiratory volume variations, gastric oscillations, vasomotion (Tong et al., 2019), and vascular structure (Aso et al., 2017) including specifically deoxyhemoglobin density (Aso et al., 2019).

With growing evidence that hemodynamic timing is a useful clinical metric in pathologies such as moyamoya disease (Donahue et al., 2016; Jahanian et al., 2018), carotid stenosis (Chang et al., 2013), stroke (Ni et al., 2017; Siegel et al., 2016), glioma (Cai et al., 2022), and Alzheimer’s (Holmes et al., 2020), we must be cautious when interpreting fMRI timing information that is derived from different stimuli or conditions. In some experimental circumstances, it may be desirable to deliver a controlled and substantial vasodilatory stimulus to enhance the effect size of the corresponding BOLD signal changes and make characterization of hemodynamic timings more robust. Researchers use complex gas inhalation challenges (Liu et al., 2019) as well as simple breathing tasks (e.g., breath-holding) in order to modulate blood CO2 levels and evoke a large systemic vascular response. Breathing-tasks have been used successfully to map the variability of BOLD hemodynamic timings in both healthy and patient populations (Bright et al., 2009; Chang et al., 2008; Geranmayeh et al., 2015; Magon et al., 2009; Moia et al., 2021; Pinto et al., 2016; Raut et al., 2016; Sousa et al., 2014). We have previously demonstrated that the timing of the BOLD CVR response is more robustly characterized when a large systemic vascular response is evoked by a breathing challenge, compared to when only spontaneous fluctuations in arterial CO2 levels are present (Bright and Murphy, 2017; Stickland et al., 2021). Still, researchers have had success with mapping BOLD hemodynamic timings in healthy and patient populations using resting-state fMRI data in the absence of evoked responses (Cai et al., 2022; Hu et al., 2021; Liu et al., 2017; Siegel et al., 2016; Tong et al., 2019). Despite the anticipated lower signal-to-noise ratio of the vascular effects, this resting-state experimental design is simpler and therefore suitable for a wider range of populations and clinical or research contexts.

However, we anticipate that the mechanisms underpinning voxel-wise and global hemodynamics in these two conditions, breathing-task and resting-state fMRI, may differ, making it unclear whether hemodynamic timings derived from these two types of data provide interchangeable physiological information. Although it is reasonable to assume that the majority of information present in the fMRI time-series during a breath-holding task reflects the associated evoked CVR response, there will be some neural activity associated with the voluntary performance of the task, in addition to a baseline level of intrinsic neuronal fluctuations. Similarly, intrinsic low-frequency fluctuations in resting-state fMRI signals are frequently used to characterize neural activation patterns of intrinsic brain networks, but it is well established that many physiologic phenomena, both neural and vascular, influence blood flow and BOLD fMRI signals in uncontrolled ways (Caballero-Gaudes and Reynolds, 2017; Liu, 2013, 2016; Murphy et al., 2013; Tong et al., 2019). Finally, systemic physiological drivers of LFOs other than CO2 will impact both conditions (Aso et al., 2019). Thus, these two types of fMRI data – breathing-task and resting-state – contain diverse types of neural and vascular information, albeit with different weightings.

In this paper, we estimate relative hemodynamic timings across the brain in BOLD fMRI data from healthy participants. We define a voxel’s hemodynamic “delay” as the time a BOLD amplitude change occurs with respect to an amplitude change in a reference signal (also called a ‘probe regressor’). These are relative hemodynamic delays, representing the timing response of one brain region relative to another; this is therefore not directly capturing arterial transit times, though regional differences in transit times may affect these estimates. Rather than emphasizing local neural activity patterns, which may or may not relate to all brain regions, the probe regressor is commonly a physiological signal reflecting a more systemic fluctuation. The probe may be derived from an external recording related to cardiac or respiratory-related phenomena. Alternatively, a data-driven probe regressor can be used, typically the average BOLD signal in brain or gray matter tissue filtered to a certain frequency band, which is thought to mostly reflect these global systemic drivers (Aso et al., 2017; Liu et al., 2017; Tong et al., 2019). Previous work has also shown maps of hemodynamic delays to have good inter-session reproducibility and be more stable across different experimental designs when using a recursive tracking method for modeling (as implemented in this paper), versus a simple seed-based cross correlation approach (Aso et al., 2017).

Using this data-driven strategy, we assess the agreement of voxel-wise hemodynamic delay estimates derived from data that contained a breath-hold task and from resting-state data, to understand how similar these delay estimates are within an individual. To further interpret our assessment of the agreement, we applied amplitude thresholds (Chang et al., 2008; Fesharaki et al., 2021) and different levels of spatial smoothing of the input data, hypothesizing that signal-to-noise factors would modulate the absolute agreement of hemodynamic delay estimates. We also investigated which areas of the brain (if any) had the highest agreement, to help provide anatomical insight into what may be driving agreement or disagreement. Comparing these two paradigms may shed light on their different physiological drivers, or suggest they provide similar information under certain methodological constraints.

2. Methods

The data from this study unfortunately cannot be made openly available due to restrictions of the ethical approval that they were collected under. The visual instructions for the breathing tasks, displayed during scanning, were created with PsychoPy code (Stickland et al., 2022) and the hemodynamic timing estimates were produced with the Rapidtide Toolbox v2.0.9 (Frederick, 2016–2022). Specific details on how these coding repositories were used, as well as other primary analysis code for this project, have been organized into this GitHub repository: github.com/BrightLab-ANVIL/Gong_2023.

2.1. Data acquisition

This study was reviewed by Northwestern University’s Institutional Review Board. Written informed consent was obtained from all subjects. Nine healthy volunteers were included in the study (6 females, mean age = 26.2 ± 4.1 years).

A 3T Siemens Prisma MRI system with a 64-channel head coil was used to collect neuroimaging data. The functional T2* -weighted acquisitions run during the hybrid breathing-task and resting-state protocol (Fig. 1A) were gradient-echo planar sequences provided by the Center for Magnetic Resonance Research (CMRR, Minnesota) with the following parameters: TR/TE = 1200/34.4 ms, FA = 62°, Multi-Band (MB) acceleration factor = 4, 60 axial slices with an ascending interleaved order, 2 mm isotropic voxels, FOV = 208 × 208 mm 2, Phase Encoding = AP, phase partial Fourier = 7/8, Bandwidth = 2290 Hz/Px. One single-band reference (SBRef) volume was acquired before the functional T2* -weighted acquisition (the same scan acquisition parameters without the MB acceleration factor) to facilitate functional realignment and masking. A whole brain T1-weighted EPI-navigated multi-echo MPRAGE scan was acquired, adapted from (Tisdall et al., 2016), with these parameters: 1 mm isotropic resolution, 176 sagittal slices, TR/TE1/TE2/TE3 = 2170/1.69/3.55/5.41 ms, TI = 1160 ms, FA = 7°, FOV = 256 × 256, Bandwidth = 650 Hz, acquisition time of 5 min 12 s, including 24 reacquisition TRs. The three echo images were combined using root-mean-square.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0001.jpg

Panel A: The BOLD time-series averaged across gray matter (GM) for each subject (thinner blue lines) and the group average (thicker line). This time-series data was cut into two segments: one segment, BH+REST, includes a breathing task (solid blue bar on x-axis) followed by REST; the second segment is only REST. 390 vol amount to 7 min 48 s (TR=1.2 s). Panel B: The output maps, from Rapidtide, for one example subject. The top row shows the correlation amplitude at each voxel (maximum correlation coefficient between each voxel time-series and probe regressor). The bottom row shows the relative delays at each voxel (time in seconds when the maximum correlation occurred). Panel C: Spatial correlation, in GM voxels, between hemodynamic delay times (panel B, bottom) from the BH+REST and REST data segments of an example subject. The correlation amplitude (panel B, top) is used to threshold the delays; the figure legend shows 3 examples of amplitude thresholding we applied. Each dot represents a voxel passing the threshold for both data segments, e.g., 0 = hemodynamic delay for a specific voxel must be in GM and have a correlation amplitude greater than 0 to be included. The lines are linear regression lines. For comparison, the y = x line is also shown.

Previous work with similar tasks (Bright et al., 2009; Lipp et al., 2015) was used as guidance for the breath-hold (BH) task timings. Subjects practiced the BH tasks outside the scanner with the researcher (R.C.S) prior to scanning. During the scan, subjects performed three BH cycles (paced breathing and a 15-second end-expiration hold) followed by an 8-minute period of rest with visual fixation. Two equal-length data segments (Fig. 1A) were analyzed: one comprised the breathing task followed by rest (BH+REST), whereas the other just included rest (REST). These neuroimaging data were used within previously published work (Stickland et al., 2021).

2.2. Data analysis

2.2.1. Data preprocessing

A custom shell script grouped AFNI (Cox, 1996) and FSL (Jenkinson et al., 2012; Li et al., 2016; Smith et al., 2004; Woolrich et al., 2009) commands. DICOMS were converted to NIFTI format with dcm2niix (Li et al., 2016). For motion correction of the fMRI dataset, AFNI’s 3dvolreg was performed with the SBRef as the reference volume; six motion parameters (three translations, three rotations) were extracted, demeaned, and saved. The first 10 vol were discarded to allow the signal to achieve a steady state of magnetization and then brain extraction (Smith, 2002) was performed using FSL. The T1-weighted dataset was processed with the fsl_anat function, involving brain extraction, bias field correction, and tissue segmentation with FAST (Zhang et al., 2001). A GM tissue mask was subsequently created by thresholding the partial volume estimate image to 0.5 and transforming it to the subject-specific fMRI space with FSL’s FLIRT command (Jenkinson et al., 2002; Jenkinson and Smith, 2001). The fMRI dataset was portioned into two segments (BH+REST and REST), each containing 390 vol (Fig. 1A). Finally, detrending and denoising were performed via linear regression using AFNI’s 3dTproject command, including the removal of Legendre polynomials up to the 4th degree to remove drifts in the signal and the six motion parameters. This single-subject preprocessed fMRI dataset, divided into two equal-length segments, was the input to the Rapidtide function, explained below.

The MNI structural atlas (Collins et al., 1995; Kötter et al., 2001), which parcellates the brain into 9 anatomical regions (caudate, cerebellum, frontal lobe, insula, occipital lobe, parietal lobe, putamen, temporal lobe, thalamus), was linearly transformed to fMRI subject space using FSL FLIRT, then masked to the subject’s GM. This transformation was achieved using nearest neighbor interpolation and 12 degrees of freedom and concatenating previously created registration matrices (the functional to structural space registration matrix and the structural to MNI registration matrix were previously created from running the fsl_anat function).

2.2.2. Estimation of relative hemodynamic timings

The Rapidtide toolbox (v2.0.9) a suite of Python programs, was used to create maps of relative hemodynamic timing. Rapidtide implements an algorithm called RIPTiDe (Regressor Interpolation at Progressive Time Delays) that performs rapid time delay analysis in filtered fMRI data to find lagged correlations between each voxel’s time-series and a reference time-series, referred to as a “probe regressor” (Frederick, 2016–2022). We used a data-driven probe regressor, initially defined as the mean GM time-series. We chose to use the Rapidtide toolbox to estimate these relative hemodynamic timings because: (i) it is a citable open-source software package, (ii) it has been demonstrated to be successful when applied to resting-state data (Champagne et al., 2022; Erdogan et al., 2016) and (iii) with a data-driven probe regressor, the resting-state and breath-hold data can be analyzed with an identical pipeline.

The Rapidtide commands are included in the GitHub repository listed earlier in this manuscript, and the input arguments used are summarized in Supplementary Table 1. These commands were run separately on the BH+REST and REST fMRI data segments for all subjects. Two Rapidtide output maps were extracted for further analysis (Fig. 1B): the hemodynamic delay (temporal offset of the maximum correlation with the probe regressor, in seconds) and the correlation amplitude (maximum similarity function value, Pearson’s r). Rapidtide resamples the probe regressor in each refinement step and produces correlation estimates on a continuous scale within the search range given. We chose a search range of −15 to +15 s to comfortably cover expected physiological ranges in a healthy cohort alongside using a data-driven probe regressor with no large offset. In our further analyses of delay values, we did not include any values at the boundary of this search range, deeming these not a true optimization and likely not physiologically plausible.

2.2.3. Hemodynamic delay comparison

The voxel-wise agreement between delay values from BH+REST and REST was calculated for all GM voxels using linear regression analysis, via the polyfit function in MATLAB (MATLAB and Statistics Tool-box Release 2021b, The MathWorks, Inc., Natick, Massachusetts, United States). This was repeated for subsets of GM voxels with increasing correlation amplitude thresholds applied (step = 0.1, max 0.7): a voxel’s correlation amplitude value had to pass the threshold in both BH+REST and REST data segments of a given subject to have the delay values included.

Three metrics were summarized: the spatial correlation coefficient (representing model fit), the slope of the linear regression model (representing absolute agreement in values), and the percentage of voxels remaining after amplitude thresholding with respect to the entire GM. When the amplitude threshold is labeled as ‘None’ (e.g., x-axes in Fig. 2), the percentage of voxels remaining is not at 100% because Rapidtide fails to produce a hemodynamic delay estimate for select GM voxels, resulting in the initial exclusion of these voxels from analysis. Fig. 1C demonstrates the relationship between delay values at three different amplitude thresholds, for one example subject. All three metrics were then summarized across subjects using boxplots. The three metrics were also investigated within subregions of GM voxels. Using the MNI atlas described above, (Section 2.2.1), the same linear regression was performed for each GM atlas region (or parcel), then a median value across subjects was calculated for each metric, in order to easily visualize and compare the pattern of results between GM subregions.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0002.jpg

Delay agreement between BH+REST and REST data segments within GM voxels, for each amplitude threshold. The boxplots summarize across subjects. The left plot shows the percentage of GM voxels satisfying each amplitude threshold in both data segments, where ‘None’ refers to no amplitude thresholding applied, and the subsequent numbers (0 to 0.7) refer to the amplitude threshold the voxel must exceed in order to be included in the assessment of delay agreement. The middle and right plots summarize the delay agreement for each amplitude threshold, showing the spatial correlation coefficients (middle plot) and slopes (right plot).

2.2.4. Impact of spatial smoothing

Prior research using a lagged-GLM approach to map regional hemodynamic timings showed that these results were more robustly achieved in breath-hold versus resting-state data (Bright and Murphy, 2017; Stickland et al., 2021), likely due to the reduced signal-to-noise of respiratory-evoked vascular responses in resting-state data. To evaluate the role of signal-to-noise (SNR) limitations in the current results, we enhanced the degree of spatial filtering applied to the data prior to hemodynamic delay estimation, from σ = 1 mm (FWHM 2.35 mm, Rapidtide recommended setting for smoothing) to σ = 2.13 mm (5 mm FWHM). The same three metrics were calculated and compared with the original less-smoothed results.

2.2.5. Impact of overlapping time-series data

In the above analysis, the breathing task data segment and the resting-state data segment were extracted from the same functional scan acquisition for each of the subjects. In order to maximize and match the degrees of freedom between data segments, there is substantial overlap in the time-series information between the BH+REST and REST data segments (as indicated in Fig. 1A). It is unclear whether this time-series overlap is critical to our observations of delay mapping. Therefore, in order to determine the effect of time-series overlap, an additional comparison of the BH+REST data segment against the REST data segment from another scan was conducted (Specifically, we isolated the resting state segment from the CDB+REST scan acquired in the same session, where “cued deep breathing” (CDB) is an alternative breathing task to BH, described in greater detail in Section 4. 3.)

3. Results

3.1. Hemodynamic delay agreement between breathing-task and resting-state data

Fig. 1B displays amplitude and delay maps for one example subject, and maps for all subjects are presented in Supplementary Figure 1. For both BH+REST and REST data segments, correlation amplitude values were generally positive and show tissue contrast between GM and white matter (WM). Voxels with negative correlations, or where Rapidtide was not able to output a value, appear predominantly in or adjacent to the cerebrospinal fluid (CSF) and sometimes in WM; this is more notable for the REST data segment. The REST data segments also display lower correlation amplitudes across the brain compared to BH+REST, as anticipated. Delay maps show good GM/WM contrast for both data segments, although some voxel clusters of extreme values are present (see Discussion).

Fig. 2 shows the summary of the three metrics used to assess hemodynamic delay agreement at different correlation amplitude thresholds: percentage voxels remaining, spatial correlation, and linear slope. As the amplitude threshold increases, fewer voxels remain, and the spatial correlation and slope of agreement relating the BH+REST vs REST hemodynamic delays both increase in a non-linear manner (Fig. 1C, Fig. 2). In other words, the agreement between GM hemodynamic delay times in REST and BH+REST increases with stricter amplitude thresholding, as shown by the higher spatial correlation and the slope approaching a value of 1. Scatterplots comparing BH+REST and REST delay values for select amplitude thresholds are provided for all subjects (Supplementary Figure 2).

The parcellation analysis showed similar trends in each brain region examined (Fig. 3). For each brain region, the percentage of voxels remaining (for each data segment separately and their union), the spatial correlation and the slope all changed similarly with increasing amplitude threshold.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0003.jpg

Delay agreement between BH+REST and REST data segments within GM parcel voxels, for each amplitude threshold. The three plots on the top row are displayed in the same format as Fig. 2 but display the median value across subjects for each GM parcel (colored lines) alongside the median value for total GM (black line), for visualization simplicity. Additionally, the bottom row shows the percentage of voxels remaining for both data segments separately. The highest amplitude threshold that is summarized is capped at 0.5 because, for some GM parcels and some subjects, too few voxels pass the higher amplitude thresholds.

After visualizing the spatial location of voxels passing high amplitude thresholds (0.6, 0.7), we noted their proximity to the sagittal sinus and other regions that are expected to have large venous blood vessel contributions to the fMRI signal (Fig. 4). Based on this observation, we hypothesized that higher amplitude thresholds may isolate voxels reflecting larger veins rather than tissue. The voxel-wise Mean Signal Intensity (MSI) is partially reflective of local venous weighting: venous blood induces large signal dephasing, resulting in lower MSI across the functional BOLD-weighted images. These voxels may drive agreement between BH+REST and REST delay values, without faithfully representing tissue physiology.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0004.jpg

Voxels passing the 0.6 amplitude threshold in both BH+REST and REST data segments overlaid on mean signal intensity (MSI) maps, and histograms showing MSI distributions for different amplitude threshold bands, for two example subjects and group-level summary. Example subject maps are in subject space and the group result maps are in MNI space.

To address this, we investigated the relationship between MSI and correlation amplitude threshold, hypothesizing the voxels passing the higher thresholds in both BH+REST and REST data segments would have lower MSI than the other thresholds, potentially due to the presence of larger veins. To visualize, MSI maps were overlayed with voxels passing the 0.6 amplitude threshold for both single subject and group data (Fig. 4). Using the ksdensity function in MATLAB, the histograms of MSI values of voxels corresponding to each of the amplitude threshold bands 0–0.2, 0.2–0.4, 0.4–0.6, and larger or equal to 0.6 were extracted (Fig. 4, with results from all subjects shown in Supplementary Figure 3). We observed that voxels with the highest correlation amplitudes in both data segments (r > 0.6) did generally reflect lower MSI values than other GM voxels, supporting our hypothesis that these voxels, at least in part, contain larger venous vessels.

To determine whether these voxels with very high correlation amplitudes drive the pattern of results presented in Fig. 2, the original analysis was repeated for correlation amplitude bands rather than lower thresholds: 0 to 0.2, 0.2 to 0.4, and 0.4 to 0.6. However, with these amplitude threshold bands, only voxels within the threshold band for both data segments can be included, which ultimately resulted in only ~25% of total GM voxels being compared in subsequent analyses (Supplementary Figure 4). This could have originated from the higher SNR of the BH+REST data segment compared with that of the REST task, such that global vasodilatory effects led to systematically higher correlation amplitude values. We therefore standardized correlation amplitudes in each data segment into percentile scores, and examined voxels that were in the same quartile for both data segments in an attempt to include more GM voxels in the subsequent analysis. Scatterplots showing these relationships for all subjects are provided in Supplementary Figure 5. Figure 5 shows the summary of the three metrics used to assess hemodynamic delay, revealing similar trends to the original results: the spatial correlation increases and the slope increases, towards 1, with each increasing amplitude percentile threshold band. Note that the percentages of voxels included in these threshold bands are still small after matching by percentile scores, indicating that many voxels show very different correlation amplitudes in the two data segments.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0005.jpg

Delay agreement between BH+REST and REST data segments within GM voxels, for each amplitude percentile threshold band. The boxplots summarize across subjects. These results are displayed in the same format as Fig. 2 (left: percentage of GM voxels satisfying each amplitude threshold band, middle: spatial correlation coefficients, right: slopes) but instead of filtering the voxels with a single lower amplitude threshold, amplitude percentile threshold bands are used (a voxel must be within the percentile threshold band in both data segments to be included).

3.2. Impact of spatial smoothing

Under additional spatial smoothing, we observe a widespread increase in correlation coefficient amplitudes (Fig. 6) in both the BH+REST and REST data, which most likely reflects the improved signal-to-noise ratios achieved through this preprocessing step. Maps generated from the smoothed data from all subjects are presented in Supplementary Figure 6, and scatterplots comparing the resulting BH+REST and REST delay values for select amplitude thresholds and threshold bands are provided for all subjects (Supplementary Figures 7 and 8). Fig. 7A demonstrates the effects of spatial smoothing on agreement between relative hemodynamic delay estimates, with the original results (median values from Fig. 2) provided for comparison. Due to the increased correlation coefficient amplitudes in both data segments, the percentage of voxels remaining increased after additional smoothing. However, the spatial correlations of hemodynamic delays observed at each threshold are mostly unchanged. There is trend for better agreement (slope closer to 1) at lower amplitude thresholds in the additionally smoothed data (Fig. 7A), however this is not readily interpretable given the inter-subject variability and the small sample size available. For the amplitude threshold band analysis after additional smoothing (Fig. 7B, Supplementary Table 2), we see a different pattern in the voxel percentages in each threshold band. This corresponds to our observation that, with additional smoothing, more voxels pass higher amplitude thresholding. For every amplitude threshold band, we observe an increase in group median correlation amplitude and slope values after additional smoothing.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0006.jpg

Parameter maps generated with Rapidtide, for one example subject (sub01), before and after changing the level of smoothing. The top row displays the relative delays at each voxel (time in seconds when the maximum correlation occurred). The bottom displays the correlation amplitude at each voxel (maximum correlation coefficient between each voxel time-series and probe regressor).

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0007.jpg

Panel A and B summarize amplitude threshold and amplitude percentile threshold band results, respectively, comparing the results with increased smoothing (FWHM 5 mm) to the original results which had less smoothing (FWHM 2.35 mm). The original results are shown as a black line for reference and duplicate the group median results from Figs. 2 and and55.

3.3. Impact of overlapping time-series data

The BH+REST and REST data segments are derived from the same scan and contain substantial overlapping time-series information. We examined the impact of this overlap on hemodynamic timing agreement, by comparing the BH+REST data segment with the REST segment from an entirely distinct functional acquisition (described in more detail in Section 4.3). The results are presented in Supplementary Figure 9. We observe the same effects as in our original analyses, suggesting that the agreement in hemodynamic timings is not entirely driven by the overlap in data segments.

4. Discussion

In BOLD-fMRI data containing a breath-hold task, we assume the globally driven vasodilatory effect dominates local hemodynamics, leading to a stronger relationship between a GM “probe” time-series and each voxel time-series. In comparison, in data containing only spontaneous resting-state fluctuations, the signal amplitude and timing characteristics at each voxel will be more region-specific, influenced more by neurovascular factors which are not globally driven. In this study, we modeled resting-state and breathing-task data segments with an identical analytical pipeline, to assess whether they produced equivalent estimates of relative hemodynamic timings across cortical brain regions. As we applied stricter correlation amplitude thresholds (correlation between the voxel time-series and the probe regressor), including a smaller percentage of the total GM voxels with each new threshold applied, the spatial agreement of hemodynamic delays between BH+REST and REST increased. Broadly, the same pattern of results was observed both when assessing all GM together and specific GM regions separately. We further observed that voxels passing the highest amplitude thresholds were more likely to be in closer proximity to venous structures; when removing these voxels from our analysis, we found similar trends to the original results, suggesting the delay agreement was not only driven by voxels proximal to large veins. Finally, applying additional smoothing to the input data boosted SNR and enhanced common variance shared across neighboring voxels; this produced widespread increases in correlation coefficient amplitudes, allowing more voxels to be included in the assessment of delay agreement. After additional smoothing, the pattern of agreement in delay values (as stricter thresholds were applied) did not change, yet some improvement in delay value agreement was observed at the group average level. The subtleties of these results, alongside their limitations and implications, are further discussed.

4.1. Estimating hemodynamic delays

Our implementation of the Rapidtide toolbox estimated a delay value for the majority of voxels in the brain, with the only inputs being the fMRI data and a GM tissue mask generated from a T1-weighted image. Therefore, this is a feasible and accessible way for researchers to map hemodynamic delays and can be implemented with the imaging data alone (i.e., without the need for additional equipment or recordings). The majority of the Rapidtide settings we applied are default or recommended by the toolbox developers with some input arguments adapted specifically for the nature of our dataset and questions, as described in the methods section and Supplementary Table 1. In Tong et al., 2019, the Rapidtide authors explain some of the features of using a global mean signal as the probe regressor: it is simple to implement, does not require an externally derived regressor, and resembles the signal fluctuations of the superior sagittal sinus (aligned with our observations in Fig. 4). This will be true for the GM average signal as well, which we used to define and refine the probe regressor in this study. We chose to focus on GM to match the voxels where agreement of hemodynamic delay values was assessed, acknowledging that GM tissue has the most relevance to fMRI analyses that rely on neurovascular coupling. Furthermore, we did not want timing differences between GM and WM voxels to dominate our assessment of agreement between data segments.

In fMRI analyses, the global signal has traditionally been used within the context of denoising, with the reasoning that it captures global “noise” rather than local neural “signal”. However, its use is controversial considering the multitude of sources that can contribute to this signal (Liu et al., 2017). With the same logic that removing the global signal allows for the remaining variance to more closely resemble neural-driven mechanisms, specifically studying the global signal and its relationship to local signals can tell us more about non-neural hemodynamics. Tong et al. (2019) discuss the challenges of using a data-driven mean signal regressor for both purposes; it is temporally “blurred” because it contains contributions from voxels with a range of delays and will also be influenced by local neural activity. To address these challenges Rapidtide implements a cross-correlation approach, often referred to as dynamic global signal regression in a denoising context (Erdogan et al., 2016). Within this cross-correlational approach, the Rapidtide algorithm goes through multiple analytical “passes”, a bootstrapping method to improve the probe regressor and make it more closely resemble non-neural processes (see the flowchart in Fig. 5 of Tong et al., 2019). Further features, such as temporal filtering and iterative depseckling (using information from neighboring voxels to correct potentially bad delay estimates) are also applied. Please refer to Tong et al. (2019) and the Rapidtide documentation (Frederick, 2016–2022) for more details and a more in-depth discussion of how voxel-wise hemodynamic delays are estimated and interpreted within this software package.

From visual inspection, the Rapidtide amplitude and delay maps broadly agree with previous CVR literature in healthy cohorts using resting-state or breathing-task designs, showing contrast between GM and WM tissue (Moia et al., 2021; Stickland et al., 2021; Tong et al., 2019), with negative CVR-related correlation amplitudes scattered more commonly across WM than GM but mostly clustered adjacent to CSF (Bright et al., 2014; Thomas et al., 2013). For some subjects, the delay maps contain clusters of extreme delay values (Supplementary Figures 1 and 6). Obvious explanations are unclear, but the Rapidtide despeckling procedure, and the smoothing we applied, could have accentuated these extreme clusters. They are present in both data segments but are more prominent in the REST data segment. Considering some of these regionally specific extreme delays, we also investigated the relationship between delay values at the level of GM cortical parcels. This allowed us to assess whether it was valid to summarize the patterns of agreement in delay at the level of the whole GM, or if different brain regions displayed notably different patterns of agreement across amplitude thresholds. However, we did not observe a meaningful difference in hemodynamic delay agreement across cortical regions.

4.2. Assessing the agreement of hemodynamic delays at different amplitude thresholds

For all subjects, we see positive correlations and slopes between BH+REST and REST data segments, at every amplitude threshold. The median slope across subjects, for every amplitude threshold, is below 1. This suggests a larger range of delay estimates in the data segment containing the breathing-task compared to the data segment containing resting fluctuations only. The results and relationships we observe across all GM voxels are consistent with those observed in the smaller GM subregions. It is unclear why the caudate stands out, with a slightly different pattern compared to other regions, though it still follows the same overall trend of increasing correlation amplitude and slope across increasing thresholding. To better visualize the agreement (or disagreement) between data segments, Fig. 8 presents the ratio of the average (across subjects) correlation amplitude maps and average difference of the hemodynamic delay maps, using a visualization methodology similar to previous literature (Tong et al., 2011). It is clear that correlation amplitudes are systematically higher (ratio > 1) in the BH+REST data segment, as expected and visible in earlier figures. Interestingly, the hemodynamic delay differences are centered near zero, suggesting no systematic offset between the data segments, and the average distribution of this “error” in hemodynamic delay falls primarily within ±1TR for both GM and white matter (WM) voxels.

An external file that holds a picture, illustration, etc.
Object name is nihms-1900796-f0008.jpg

Top row: Ratio of the average correlation amplitude maps for BH+REST and REST data segments (Median: GM 1.27, WM 1.24); Bottom row: Average hemodynamic timing difference maps, BH+REST-REST (Median: GM −0.05, WM 0.04). In both rows, kernel density plots are provided to summarize gray matter (GM) and white matter (WM) values.

Further work is needed to concretely establish the reasons for (dis)agreement between hemodynamic delays measured in the breathing-task and resting-state data: are the relationships we observe driven by data quality, methodology, physiology, or a mixture of these factors? We attempted to explore these questions by assessing how the agreement changed when: (1) restricting the voxels included (via amplitude thresholds), (2) investigating specific voxels that may have a strong influence on the overall agreement and (3) changing the signal-to-noise ratio of the input data via increased spatial smoothing.

Firstly, if we examine the agreement of delay values with no amplitude threshold applied, therefore including most of the GM voxels within the linear model, we see very low correlations and slopes of 0.278 and 0.440 (median across subjects, less-smoothed and more-smoothed analysis, respectively). Moderate to strong correlations, with slopes closer to one, are only achieved when including a very small percentage of GM voxels, the voxels with the strongest positive correlation with the probe regressor. As a complementary representation of these results, we have also provided a summary of the root mean square error (RMSE) of the linear fits (Supplementary Figure 10). We show that RMSE decreases, meaning that “goodness of fit” increases, at higher amplitude thresholds, mirroring the correlation coefficient results. It is more difficult to interpret slopes when the correlation values are low, however the monotonic increase of correlation and slope, as amplitude threshold increases, is interesting to reflect on. With these results, we cannot conclude that the delay values estimated with resting-state and with breathing-task data have the same physiological meaning in every voxel; estimates may only converge when the weighting of neural, vascular and other physiological factors are the same. The improved delay agreement with higher amplitude thresholding may also reflect improved SNR, where higher correlation amplitudes make the measurement of hemodynamic delay more robust in the presence of other ongoing signal processes (“noise”).

Secondly, our results suggest that voxels with the strongest relationship with the probe regressor (in both data segments) are located close to the sagittal sinus and other regions that are expected to have large venous blood contributions. A similar observation was noted by a different group of researchers (Tong et al., 2011); they found, as have others, that the largest BOLD CVR-related amplitude changes are in large blood vessels, and that signal amplitude changes were most similar between a resting-state and a breath-hold dataset in these large vessels compared to in other regions. We went on to exclude these voxels in an attempt to assess hemodynamic delays relevant to local physiology, not just global vascular drainage. Without the inclusion of these high-amplitude voxels (within the upper quartile), we still observe a monotonic increase of spatial correlation and slope with each increasing amplitude percentile threshold band. This provides further evidence indicating that voxels with large venous vessels drive some of the agreement seen, but likely not all. Larger samples, and data acquisitions more tailored towards identifying venous weighting, would be needed to understand if this is meaningful.

Finally, increased smoothing of the input data increased the amplitude correlation coefficients of GM voxels on average (i.e., increased their correlation strength with the probe regressor). This is likely because of a reduced contribution from neurovascular fluctuations specific to each voxel, as well as other sources of local fMRI noise, leaving more shared features between neighboring voxels. This increased smoothing allowed more GM voxels to be included in the assessment of delay agreement for every amplitude threshold and manifested in similar or slightly improved hemodynamic delay agreement between the resting-state and breathing-task datasets. Therefore, this suggests that SNR factors do play a role in the agreement of delay estimates, and this is also consistent with the idea that higher agreement is observed in voxels with a stronger relationship to the global probe regressor. However, smoothing the data also amplified some unusual disagreement in hemodynamic delay observed in select subjects. For example, the scatterplots for Subject 3 in Supplementary Figures 2A and 7 demonstrate that there are select groups of voxels with near-zero hemodynamic delay in the BH+REST data but much more extreme relative delays in the REST data segment. (Note that this phenomenon appears limited to a small number of voxels, and it is counter to our observation, described above, that breathing-task data generally produces a wider range of delay values across all GM.) Such disagreement could suggest that the Rapidtide algorithm is problematically impacted by other ongoing signal fluctuations present in resting-state data that have distinct dynamics to low-frequency systemic vasodilatory effects and are coordinated across multiple contiguous voxels. In contrast, the higher SNR and greater specificity to CO2-related vasodilatory effects that are achieved by having a breathing task may prevent this fitting challenge, irrespective of spatial smoothing.

4.3. Alternative breathing tasks

As an alternative to the BH task which induces hypercapnia, the neuroimaging dataset also included a cued deep breathing (CDB) task in the same participants, which induces hypocapnia via hyperventilation, described in (Stickland et al., 2021). Two repetitions of fast deep breaths were cued, and an 8-minute rest period followed this CDB task. As with the BH+REST dataset, two equal-length data segments were extracted. We replicated the main analysis with this alternative breathing task (the agreement of hemodynamic timings between CDB+REST and REST, across amplitude thresholds), shown in Supplementary Figures 2B and 11. Except for one subject being removed from group summaries as an outlier, we note the same trends and key findings reported in the main manuscript with the BH+REST dataset. Although our observations may be fairly robust to the breathing task used, we note that this may not be true in patients with vascular pathology; altered baseline vascular tone can differentially impact hypercapnic and hypocapnic cerebrovascular reactivity, including both the magnitude and delay of the vascular response (Bright et al., 2011; Zhao et al., 2009).

4.4. Alternative reference signals

The complexity of the vascular tree, and the heterogeneous venous-weighting of BOLD fMRI contrast, makes it difficult to disentangle the fundamental biophysical meaning of the hemodynamic delays that we are measuring in our analysis. The reference time-series (e.g., probe regressor) used in our analysis is recursively derived from the GM mean time-series, which will reflect a mixed sampling of vascular sources. In the existing literature using similar delay mapping methodology, other groups have used a more purely venous reference time-series, derived from the superior sagittal sinus (SSS) (Aso et al., 2017; Tong et al., 2016; van Niftrik et al., 2016). Imitating their methods, we manually delineated an ROI in the SSS of each subject’s T1-weighted anatomical scan, registered this ROI mask to their functional data, and used this mask to define a probe regressor in Rapidtide. The correlation amplitude maps showed strong correlations and high agreement between GM and SSS results when compared within a GM ROI, particularly for the BH+REST segment (Supplementary Table 3). For delay map agreement, there was more subject variability, but on average these correlations are still moderate to strong, and again stronger for the BH+REST data segment. To bring context to strength of these correlations – no amplitude thresholding was applied within the GM ROI for the results presented in Supplementary Table 3, therefore this suggests greater similarity for maps produced with different reference signals (but the same data segment) compared to maps produced with the same reference signal (but different data segments), as shown in manuscript Fig. 2 for GM with no amplitude threshold applied. We then replicated the original analyses described above, reproducing Fig. 2 using the new SSS reference time-series results (Supplementary Figure 12). The trends we observe in this modified analysis (decrease in voxel percentage remaining, increase in Fisher’s Z spatial correlation, and increase in slope as the amplitude threshold is increased) are almost identical to those in the original results using a GM region to define the reference timeseries. Therefore, this change in reference time-series does not influence our main conclusions regarding hemodynamic delay timings in BH+REST versus REST data segments.

4.5. Alternative approaches to estimating hemodynamic delays

There are alternative approaches to estimating hemodynamic delays with BOLD-fMRI data, each with their pros and cons relating to complexity of acquisition and analysis. Using probe regressors external to the brain such as end-tidal CO2 or near-infrared spectroscopy (Pinto et al., 2016; Tong et al., 2019), instead of a data-driven regressor, results in the regression coefficients having a more specific physiological interpretation and being easier to compare across studies. This is a more important priority when scaling BOLD amplitude changes, rather than timing (Zvolanek et al., 2023); timing estimates are often normalized to a common spatial reference (i.e., median timing across GM voxels, as done in this paper) even when an external probe regressor is used. Furthermore, we previously found that when using an external end-tidal CO2 regressor within a lagged general linear model framework it can be challenging to get trustworthy estimates of delay in resting-state data alone (Bright et al., 2017; Stickland et al., 2021).Other modeling approaches exist to extract timing information, for example, Fourier based analyses (Pinto et al., 2016), but without a signal model for the resting-state data, this cannot be applied in the same way as with a breathing task or gas delivery design with a specific task frequency. Compared to our task paradigm, previous papers using breathing tasks will typically employ more repeated breath-hold cycles to improve the modeling of vascular reactivity mechanisms, and it is also recommended to use multi-echo BOLD fMRI data with spatial ICA-denoising to increase BOLD signal-to-noise ratio and more effectively remove motion artifacts (Cohen and Wang, 2019; Moia et al., 2021). In our data, we regressed out the six motion parameters prior to running Rapidtide, however we cannot conclude that motion is fully accounted for. These factors should be considered to improve future comparisons of BOLD timing estimates from resting-state and breathing-task designs.

4.6. Limitations and further considerations

The two data segments compared in this manuscript, BH+REST and REST, overlapped in time (Fig. 1A). It could be argued that this would lead to artificial agreement when trying to compare ‘breathing-task’ data to ‘resting-state’ data; though this may be the case (despite contrary evidence presented in Section 3.3), our results show that good agreement of hemodynamic timings is only achieved under relatively strict analytical constraints, even with this overlap. Looking at Supplementary Figures 1 and 6, there are many clusters of voxels that show divergent hemodynamic delay estimates in BH+REST compared to REST (for example see subjects 07, 09, and 10), and much of the analysis we performed attempted to understand this disagreement. The advantages of using overlapping time segments include fewer confounding factors due to the participant being in a different physiological state. Furthermore, using two different sections of the same scan means sources of error driven by volume registration are minimized, making the voxel-wise comparisons more valid. This allowed us to more confidently interpret any differences in hemodynamic timing estimates as driven by replacing a resting-state portion of the dataset with a hypercapnic breathing task portion.

With our data and analysis methodology, we cannot specifically quantify arterial arrival times or blood flow changes in response to a known change in arterial blood gas tensions. However, when accounting for voxel-wise delays in fMRI analyses, and bringing clinical insight to regional cerebrovascular dysfunction, it is likely that relative hemodynamic timings are sufficient. We also anticipate that concurrent autoregulatory phenomena are present in both data segments, and that blood pressure fluctuations may be particularly enhanced in the BH+REST data segment and time-locked to the breathing task. Future work would benefit from the addition of beat-to-beat blood pressure monitoring to better characterize these changes and facilitate the differentiation of autoregulatory processes from neurovascular coupling or cerebrovascular reactivity effects. Furthermore, the short breath-hold tasks used in this study are known to drive mild hypoxia in addition to mild hypercapnia, and hypoxia may have distinct influence on the cerebrovasculature. The two effects are profoundly correlated in breath-holding fMRI data (Bright and Murphy, 2013), and difficult to tease apart in correlative time-series analysis. Thus, in our BH+REST segment, there are several concurrent effects causing BOLD signal changes: 1) vasodilation caused by increased PaCO2, 2) vasodilation caused by reduced PaO2, and 3) reduced PaO2 itself, all caused by the apnea. The two vasodilatory responses will both incur a BOLD signal increase throughout most of cortical GM, whereas reductions in PaO2 may lead to BOLD signal decreases. We note, however, that the vasodilatory effect of mild hypoxia is anticipated to be small (Brugniaux et al., 2007). Although we do not differentiate the contributions of these three mechanisms to our observations, we can make several inferences based on our data and the literature. First, Tancredi and Hoge (Tancredi and Hoge, 2013) reported that breath-hold BOLD fMRI signal changes did not substantially differ from signal changes associated with controlled gas inhalation challenges, where hypoxia could be minimized. Second, we observed similar results when considering BH+REST and CDB+REST data, and the CDB task has potential hyperoxic, rather than hypoxic, effects. Still, to determine the precise contributions of concurrent physiological fluctuations on the hemodynamic delay structure throughout the brain, future work would benefit from more controlled physiological stimuli to independently modulate PaO2 and PaCO2, or other time-locked physiologic effects.

The band-pass filter applied by the Rapidtide algorithm (0.0009 – 0.15 Hz) should minimize the influence of the cardiac rhythm on our analysis, however aliased cardiac signals and other cardiac-mediated effects could still drive patterns in the lag structure particularly for participants with lower heart rates. Therefore, it may be worthwhile to specifically measure and model this alternative source of signal variation (to respiratory-mediated effects) in future studies. Furthermore, because the resting-state only data will have lower signal-to-noise compared to the data including the breathing task, these data may require additional noise removal to achieve more robust delay mapping. However, as the main comparison in this paper was to compare BH+REST with REST, we were hesitant to apply differing denoising strategies beyond the ones already implemented, particularly as a previous (single-echo) BOLD-fMRI study showed that delay maps lose both information and reproducibility when traditional ICA-denoising strategies are applied (Aso et al., 2017). Nevertheless, we do acknowledge that the absence of some of these denoising strategies in our resting-state data could have affected the level of delay agreement we observed. Future studies implementing multi-echo fMRI acquisition and denoising strategies could potentially help to mitigate this concern (Moia et al., 2021), while also adding insight into the biophysical mechanisms driving signal changes (Aso et al., 2019).

Finally, this study has a small sample size, and future work should investigate these questions further with larger and more diverse samples to determine generalizability. Additional data in patient cohorts with atypical regional hemodynamics would also help confirm the interpretation of our results. We provide many supplementary figures to complement the main results, in order to transparently communicate the full data variance in this healthy cohort.

5. Conclusions

We modeled BOLD fMRI data to estimate relative hemodynamic delay values across the brain, focusing on GM tissue. We investigated how these delay estimates change in the presence of a breathing-task versus in resting-state data. Hemodynamic delay information is a clinically meaningful marker of cerebrovascular function, and although breathing-task data may improve the robustness of delay mapping, resting-state data would not require the same degree of participant compliance and is therefore a very desirable approach. We found strong similarities between estimates from resting-state data and breathing-task data only when strict amplitude thresholds are applied, comparing a small percentage of GM voxels. A stronger agreement between delays estimated in these two different data types was seen in voxels that more closely resembled the global GM average, and localized (at least in part) to large venous structures. Caution must be taken when using voxel-wise delay estimates from resting-state and breathing-task data interchangeably throughout cortical GM, and further work is needed to more directly identify the local and global signal contributions that lead to their disagreement.

Supplementary Material

Acknowledgments

This work was supported by the Center for Translational Imaging at Northwestern University. The authors express special thanks to Dr. Blaise Frederick for his assistance implementing the Rapidtide software package.

Funding

This research was supported by the Eunice Kennedy Shriver National Institute of Child Health and Human Development of the National Institutes of Health under award number K12HD073945. JG was supported by multiple grants from the Northwestern University Office of Undergraduate Research, Robert R. McCormick School of Engineering and Applied Science, and Department of Biomedical Engineering. The conclusions, opinions, and other statements in this presentation are the authors’ and not necessarily those of the sponsoring institution.

Abbreviations:

BHbreath hold
BOLDblood oxygenation level dependent
CDBcued deep breathing
CSFcerebrospinal fluid
CVRcerebrovascular reactivity
FCfunctional connectivity
fMRIfunctional magnetic resonance imaging
GMgray matter
HRFhemodynamic response function
MSImean signal intensity
WMwhite matter

Footnotes

Declaration of Competing Interest

The authors declare no competing financial interests.

Credit authorship contribution statement

Jingxuan Gong: Conceptualization, Methodology, Software, Formal analysis, Data curation, Writing – original draft, Writing – review & editing, Visualization, Funding acquisition. Rachael C. Stickland: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – review & editing. Molly G. Bright: Conceptualization, Methodology, Resources, Writing – original draft, Writing – review & editing, Supervision, Project administration, Funding acquisition.

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.neuroimage.2023.120120.

Data and code availability statement

(Included verbatim in Section 2. Methods)

The data from this study unfortunately cannot be made openly available due to restrictions of the ethical approval that they were collected under. The visual instructions for the breathing tasks, displayed during scanning, were created with PsychoPy code (Stickland et al., 2022) and the hemodynamic timing estimates were produced with the Rapidtide Toolbox v2.0.9 (Frederick, 2016–2022). Specific details on how these coding repositories were used, as well as other primary analysis code for this project, have been organized into this GitHub repository: github.com/BrightLab-ANVIL/Gong_2023.

References

  • Aso T, Jiang G, Urayama SI, Fukuyama H, 2017. A resilient, non-neuronal source of the spatiotemporal lag structure detected by BOLD signal-based blood flow tracking. Front. Neurosci 11, 256. 10.3389/fnins.2017.00256. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Aso T, Urayama S, Fukuyama H, Murai T, 2019. Axial variation of deoxyhemoglobin density as a source of the low-frequency time lag structure in blood oxygenation level-dependent signals. PLoS One 14 (9), e0222787. 10.1371/journal.pone.0222787. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Bright MG, Bianciardi M, de Zwart JA, Murphy K, Duyn JH, 2014. Early anti-correlated BOLD signal changes of physiologic origin. Neuroimage 87, 287–296. 10.1016/j.neuroimage.2013.10.055. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Bright MG, Bulte DP, Jezzard P, Duyn JH, 2009. Characterization of regional heterogeneity in cerebrovascular reactivity dynamics using novel hypocapnia task and BOLD fMRI. Neuroimage 48 (1), 166–175. 10.1016/j.neuroimage.2009.05.026. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Bright MG, Donahue MJ, Duyn JH, Jezzard P, Bulte DP, 2011. The effect of basal vasodilation on hypercapnic and hypocapnic reactivity measured using magnetic resonance imaging. J. Cereb. Blood Flow Metab 31 (2), 426–438. 10.1038/jcbfm.2010.187. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Bright MG, Murphy K, 2013. Reliable quantification of BOLD fMRI cerebrovascular reactivity despite poor breath-hold performance. Neuroimage 83, 559–568. 10.1016/j.neuroimage.2013.07.007. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Bright MG, Murphy K, 2017. Cleaning up the fMRI time series: mitigating noise with advanced acquisition and correction strategies. Neuroimage 154, 1–3. 10.1016/j.neuroimage.2017.03.056 [Abstract] [CrossRef] [Google Scholar]
  • Bright MG, Tench CR, Murphy K, 2017. Potential pitfalls when denoising resting state fMRI data using nuisance regression. Neuroimage 154, 159–168. 10.1016/j.neuroimage.2016.12.027 [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Brugniaux JV, Hodges AN, Hanly PJ, Poulin MJ, 2007. Cerebrovascular responses to altitude. Respir. Physiol. Neurobiol 158 (2–3), 212–223. 10.1016/j.resp.2007.04.008. [Abstract] [CrossRef] [Google Scholar]
  • Caballero-Gaudes C, Reynolds RC, 2017. Methods for cleaning the BOLD fMRI signal. Neuroimage 154, 128–149. 10.1016/j.neuroimage.2016.12.018. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Cai S, Shi Z, Zhou S, Liang Y, Wang L, Wang K, Zhang L, 2022. Cerebrovascular dysregulation in patients with glioma assessed with time-shifted BOLD fMRI. Radiology 304 (1), 155–163. 10.1148/radiol.212192. [Abstract] [CrossRef] [Google Scholar]
  • Champagne AA, Coverdale NS, Allen MD, Tremblay JC, MacPherson REK, Pyke KE, Cook DJ, 2022. The physiological basis underlying functional connectivity differences in older adults: a multi-modal analysis of resting-state fMRI. Brain Imaging Behav. 16 (4), 1575–1591. 10.1007/s11682-021-00570-0. [Abstract] [CrossRef] [Google Scholar]
  • Chang C, Thomason ME, Glover GH, 2008. Mapping and correction of vascular hemodynamic latency in the BOLD signal. Neuroimage 43 (1), 90–102. 10.1016/j.neuroimage.2008.06.030. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Chang TY, Kuan WC, Huang KL, Chang CH, Chang YJ, Wong HF, Liu HL, 2013. Heterogeneous cerebral vasoreactivity dynamics in patients with carotid stenosis. PLoS One 8 (9), e76072. 10.1371/journal.pone.0076072. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Chen JJ, Gauthier CJ, 2021. The Role of Cerebrovascular-Reactivity Mapping in Functional MRI: calibrated fMRI and Resting-State fMRI. Front Physiol 12, 657362. 10.3389/fphys.2021.657362. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Cohen AD, Wang Y, 2019. Improving the assessment of breath-holding induced cerebral vascular reactivity using a multiband multi-echo ASL/BOLD sequence. Sci. Rep 9 (1), 5079. 10.1038/s41598-019-41199-w. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Cohen ER, Ugurbil K, Kim SG, 2002. Effect of basal conditions on the magnitude and dynamics of the blood oxygenation level-dependent fMRI response. J. Cereb. Blood Flow Metab 22 (9), 1042–1053. 10.1097/00004647-200209000-00002. [Abstract] [CrossRef] [Google Scholar]
  • Collins DL, Holmes CJ, Peters TM, Evans AC, 1995. Automatic 3-D model-based neuroanatomical segmentation. Hum. Brain Mapp 3 (3), 190–208. 10.1002/hbm.460030304. [CrossRef] [Google Scholar]
  • Cox RW, 1996. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput. Biomed. Res 29 (3), 162–173. 10.1006/cbmr.1996.0014. [Abstract] [CrossRef] [Google Scholar]
  • Donahue MJ, Strother MK, Lindsey KP, Hocke LM, Tong Y, Frederick BD, 2016. Time delay processing of hypercapnic fMRI allows quantitative parameterization of cerebrovascular reactivity and blood flow delays. J. Cereb. Blood Flow Metab 36 (10), 1767–1779. 10.1177/0271678X15608643. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Erdogan SB, Tong Y, Hocke LM, Lindsey KP, Frederick B, 2016. Correcting for blood arrival time in global mean regression enhances functional connectivity analysis of resting state fMRI-BOLD signals. Front. Hum. Neurosci 10, 311. 10.3389/fn-hum.2016.00311 [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Fesharaki NJ, Mathew AB, Mathis JR, Huddleston WE, Reuss JL, Pillai JJ, DeYoe EA, 2021. Effects of thresholding on voxel-wise correspondence of breath-hold and resting-state maps of cerebrovascular reactivity. Front. Neurosci 15, 654957. 10.3389/fnins.2021.654957. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Frederick BB, Rapidtide. [Computer Program] (2016-2022). Retrieved from https://github.com/bbfrederick/rapidtide 10.5281/zenodo.814990 [CrossRef] [Google Scholar]
  • Friston KJ, Fletcher P, Josephs O, Holmes A, Rugg MD, Turner R, 1998. Event-related fMRI: characterizing differential responses. Neuroimage 7 (1), 30–40. 10.1006/nimg.1997.0306. [Abstract] [CrossRef] [Google Scholar]
  • Geranmayeh F, Wise RJ, Leech R, Murphy K, 2015. Measuring vascular reactivity with breath-holds after stroke: a method to aid interpretation of group-level BOLD signal changes in longitudinal fMRI studies. Hum. Brain Mapp 36 (5), 1755–1771. 10.1002/hbm.22735. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Glover GH, 2011. Overview of functional magnetic resonance imaging. Neurosurg. Clin. N. Am 22 (2), 133–139. 10.1016/j.nec.2010.11.001, vii. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Handwerker DA, Ollinger JM, D’Esposito M, 2004. Variation of BOLD hemodynamic responses across subjects and brain regions and their effects on statistical analyses. Neuroimage 21 (4), 1639–1651. 10.1016/j.neuroimage.2003.11.029. [Abstract] [CrossRef] [Google Scholar]
  • Holmes KR, Tang-Wai D, Sam K, McKetton L, Poublanc J, Crawley AP, Mikulis DJ, 2020. Slowed temporal and parietal cerebrovascular response in patients with alzheimer’s disease. Can. J. Neurol. Sci 47 (3), 366–373. 10.1017/cjn.2020.30. [Abstract] [CrossRef] [Google Scholar]
  • Hu JY, Kirilina E, Nierhaus T, Ovadia-Caro S, Livne M, Villringer K, Khalil AA, 2021. A novel approach for assessing hypoperfusion in stroke using spatial independent component analysis of resting-state fMRI. Hum. Brain Mapp 42 (16), 5204–5216. 10.1002/hbm.25610. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Huneau C, Benali H, Chabriat H, 2015. Investigating human neurovascular coupling using functional neuroimaging: a critical review of dynamic models. Front. Neurosci 9, 467. 10.3389/fnins.2015.00467. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Issard C, Gervain J, 2018. Variability of the hemodynamic response in infants: influence of experimental design and stimulus complexity. Dev. Cogn. Neurosci 33, 182–193. 10.1016/j.dcn.2018.01.009. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Jahanian H, Christen T, Moseley ME, Zaharchuk G, 2018. Erroneous resting-state fMRI connectivity maps due to prolonged arterial arrival time and how to fix them. Brain Connect 8 (6), 362–370. 10.1089/brain.2018.0610. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Jenkinson M, Bannister P, Brady M, Smith S, 2002. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage 17 (2), 825–841. 10.1016/s1053-8119(02)91132-8. [Abstract] [CrossRef] [Google Scholar]
  • Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM, 2012. Fsl. Neuroimage 62 (2), 782–790. 10.1016/j.neuroimage.2011.09.015. [Abstract] [CrossRef] [Google Scholar]
  • Jenkinson M, Smith S, 2001. A global optimisation method for robust affine registration of brain images. Med. Image Anal 5 (2), 143–156. 10.1016/s1361-8415(01)00036-6. [Abstract] [CrossRef] [Google Scholar]
  • Kötter R, Mazziotta J, Toga A, Evans A, Fox P, Lancaster J, Mazoyer B, 2001. A probabilistic atlas and reference system for the human brain: international Consortium for Brain Mapping (ICBM). Philos. Trans. R. Soc. B: Biol. Sci 356 (1412), 1293–1322. 10.1098/rstb.2001.0915. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Leung J, Duffin J, Fisher JA, Kassner A, 2016. MRI-based cerebrovascular reactivity using transfer function analysis reveals temporal group differences between patients with sickle cell disease and healthy controls. Neuroimage Clin. 12, 624–630. 10.1016/j.nicl.2016.09.009. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Li X, Morgan PS, Ashburner J, Smith J, Rorden C, 2016. The first step for neuroimaging data analysis: DICOM to NIfTI conversion. J. Neurosci. Methods 264, 47–56. 10.1016/j.jneumeth.2016.03.001. [Abstract] [CrossRef] [Google Scholar]
  • Lipp I, Murphy K, Caseras X, Wise RG, 2015. Agreement and repeatability of vascular reactivity estimates based on a breath-hold task and a resting state scan. Neuroimage 113, 387–396. 10.1016/j.neuroimage.2015.03.004. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Liu P, De Vis JB, Lu H, 2019. Cerebrovascular reactivity (CVR) MRI with CO2 challenge: a technical review. Neuroimage 187, 104–115. 10.1016/j.neuroimage.2018.03.047. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Liu P, Li Y, Pinho M, Park DC, Welch BG, Lu H, 2017a. Cerebrovascular reactivity mapping without gas challenges. Neuroimage 146, 320–326. 10.1016/j.neuroimage.2016.11.054. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Liu TT, 2013. Neurovascular factors in resting-state functional MRI. Neuroimage 80, 339–348. 10.1016/j.neuroimage.2013.04.071. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Liu TT, 2016. Noise contributions to the fMRI signal: an overview. Neuroimage 143, 141–151. 10.1016/j.neuroimage.2016.09.008. [Abstract] [CrossRef] [Google Scholar]
  • Liu TT, Behzadi Y, Restom K, Uludag K, Lu K, Buracas GT, Buxton RB, 2004. Caffeine alters the temporal dynamics of the visual BOLD response. Neuroimage 23 (4), 1402–1413. 10.1016/j.neuroimage.2004.07.061. [Abstract] [CrossRef] [Google Scholar]
  • Liu TT, Nalci A, Falahpour M, 2017b. The global signal in fMRI: nuisance or Information? Neuroimage 150, 213–229. 10.1016/j.neuroimage.2017.02.036. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Magon S, Basso G, Farace P, Ricciardi GK, Beltramello A, Sbarbati A, 2009. Reproducibility of BOLD signal change induced by breath holding. Neuroimage 45 (3), 702–712. 10.1016/j.neuroimage.2008.12.059. [Abstract] [CrossRef] [Google Scholar]
  • Martindale J, Mayhew J, Berwick J, Jones M, Martin C, Johnston D, Zheng Y, 2003. The hemodynamic impulse response to a single neural event. J. Cereb. Blood Flow Metab 23 (5), 546–555. 10.1097/01.WCB.0000058871.46954.2B. [Abstract] [CrossRef] [Google Scholar]
  • Moia S, Stickland RC, Ayyagari A, Termenon M, Caballero-Gaudes C, Bright MG, 2020. Voxel-wise optimization of hemodynamic lags to improve regional CVR estimates in breath-hold fMRI. In: Proceedings of the 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC), pp. 1489–1492 [Abstract] [Google Scholar]
  • Moia S, Termenon M, Urunuela E, Chen G, Stickland RC, Bright MG, Caballero-Gaudes C, 2021. ICA-based denoising strategies in breath-hold induced cerebrovas-cular reactivity mapping with multi echo BOLD fMRI. Neuroimage 233, 117914. 10.1016/j.neuroimage.2021.117914. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Murphy K, Birn RM, Bandettini PA, 2013. Resting-state fMRI confounds and cleanup. Neuroimage 80, 349–359. 10.1016/j.neuroimage.2013.04.001. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Ni L, Li J, Li W, Zhou F, Wang F, Schwarz CG, Xu Y (2017). The value of resting-state functional MRI in subacute ischemic stroke: comparison with dynamic susceptibility contrast-enhanced perfusion MRI. Sci Rep, 7, 41586. 10.1038/srep41586 [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Penny W, Friston K, Ashburner J, Kiebel S, & Nichols T (2006). Statistical Parametric Mapping (SPM12) [MATLAB Toolbox]. London. Retrieved from www.fil.ion.ucl.ac.uk/spm [Google Scholar]
  • Pinto J, Bright MG, Bulte DP, Figueiredo P, 2020. Cerebrovascular reactivity mapping without gas challenges: a methodological guide. Front Physiol 11, 608475. 10.3389/fphys.2020.608475. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Pinto J, Jorge J, Sousa I, Vilela P, Figueiredo P, 2016. Fourier modeling of the BOLD response to a breath-hold task: optimization and reproducibility. Neuroimage 135, 223–231. 10.1016/j.neuroimage.2016.02.037. [Abstract] [CrossRef] [Google Scholar]
  • Raut RV, Nair VA, Sattin JA, Prabhakaran V, 2016. Hypercapnic evaluation of vascular reactivity in healthy aging and acute stroke via functional MRI. Neuroimage Clin. 12, 173–179. 10.1016/j.nicl.2016.06.016. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Siegel JS, Snyder AZ, Ramsey L, Shulman GL, Corbetta M, 2016. The effects of hemodynamic lag on functional connectivity and behavior after stroke. J. Cereb. Blood Flow Metab 36 (12), 2162–2176. 10.1177/0271678X15614846. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Sleight E, Stringer MS, Marshall I, Wardlaw JM, Thrippleton MJ, 2021. Cerebrovascular reactivity measurement using magnetic resonance imaging: a systematic review. Front. Physiol 12, 643468. 10.3389/fphys.2021.643468. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Smith SM, 2002. Fast robust automated brain extraction. Hum. Brain Mapp 17 (3), 143–155. 10.1002/hbm.10062. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Smith SM, Jenkinson M, Woolrich MW, Beckmann CF, Behrens TE, Johansen-Berg H, Matthews PM, 2004. Advances in functional and structural MR image analysis and implementation as FSL. Neuroimage 23 (Suppl 1), S208–S219. 10.1016/j.neuroimage.2004.07.051. [Abstract] [CrossRef] [Google Scholar]
  • Sousa I, Vilela P, Figueiredo P, 2014. Reproducibility of hypocapnic cerebrovascular reactivity measurements using BOLD fMRI in combination with a paced deep breathing task. Neuroimage 98, 31–41. 10.1016/j.neuroimage.2014.04.049. [Abstract] [CrossRef] [Google Scholar]
  • Stickland RC, Bright MG, Moia S, & Zvolanek KM, Breathing-Task-Visual-PsychoPy. [Computer Program] (2022). 10.5281/zenodo.7056657 [CrossRef] [Google Scholar]
  • Stickland RC, Zvolanek KM, Moia S, Ayyagari A, Caballero-Gaudes C, Bright MG, 2021. A practical modification to a resting state fMRI protocol for improved characterization of cerebrovascular function. Neuroimage 239, 118306. 10.1016/j.neuroimage.2021.118306. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Tancredi FB, Hoge RD, 2013. Comparison of cerebral vascular reactivity measures obtained using breath-holding and CO2 inhalation. J. Cereb. Blood Flow Metab 33 (7), 1066–1074. 10.1038/jcbfm.2013.48. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Thomas BP, Liu P, Aslan S, King KS, van Osch MJ, Lu H, 2013. Physiologic underpinnings of negative BOLD cerebrovascular reactivity in brain ventricles. Neuroimage 83, 505–512. 10.1016/j.neuroimage.2013.07.005. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Tisdall MD, Reuter M, Qureshi A, Buckner RL, Fischl B, van der Kouwe AJW, 2016. Prospective motion correction with volumetric navigators (vNavs) reduces the bias and variance in brain morphometry induced by subject motion. Neuroimage 127, 11–22. 10.1016/j.neuroimage.2015.11.054. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Tong Y, Bergethon PR, Frederick BD, 2011. An improved method for mapping cerebrovascular reserve using concurrent fMRI and near-infrared spectroscopy with Regressor Interpolation at Progressive Time Delays (RIPTiDe). Neuroimage 56 (4), 2047–2057. 10.1016/j.neuroimage.2011.03.071. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Tong Y, Hocke LM, Frederick BB, 2019. Low frequency systemic hemodynamic “noise” in resting state BOLD fMRI: characteristics, causes, implications, mitigation strategies, and applications. Front. Neurosci 13, 787. 10.3389/fnins.2019.00787. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Tong Y, Hocke LM, Lindsey KP, Erdogan SB, Vitaliano G, Caine CE, Frederick B, 2016. Systemic low-frequency oscillations in BOLD signal vary with tissue type. Front. Neurosci 10, 313. 10.3389/fnins.2016.00313. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • van Niftrik CH, Piccirelli M, Bozinov O, Pangalu A, Valavanis A, Regli L, Fierstra J, 2016. Fine tuning breath-hold-based cerebrovascular reactivity analysis models. Brain Behav. 6 (2), e00426. 10.1002/brb3.426. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • West KL, Zuppichini MD, Turner MP, Sivakolundu DK, Zhao Y, Abdelkarim D, Rypma B, 2019. BOLD hemodynamic response function changes significantly with healthy aging. Neuroimage 188, 198–207. 10.1016/j.neuroimage.2018.12.012. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Woolrich MW, Jbabdi S, Patenaude B, Chappell M, Makni S, Behrens T, Smith SM, 2009. Bayesian analysis of neuroimaging data in FSL. Neuroimage 45 (1 Suppl), S173–S186. 10.1016/j.neuroimage.2008.10.055. [Abstract] [CrossRef] [Google Scholar]
  • Zhang Y, Brady M, Smith S, 2001. Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Trans. Med. Imaging 20 (1), 45–57. 10.1109/42.906424. [Abstract] [CrossRef] [Google Scholar]
  • Zhao P, Alsop DC, Abduljalil A, Selim M, Lipsitz L, Novak P, Novak V, 2009. Vasoreactivity and peri-infarct hyperintensities in stroke. Neurology 72 (7), 643–649. 10.1212/01.wnl.0000342473.65373.80. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]
  • Zvolanek KM, Moia S, Dean JN, Stickland RC, Caballero-Gaudes C, Bright MG, 2023. Comparing end-tidal CO2, respiration volume per time (RVT), and average gray matter signal for mapping cerebrovascular reactivity amplitude and delay with breath-hold task BOLD fMRI. Neuroimage 272, 120038. 10.1016/j.neuroimage.2023.120038. [Europe PMC free article] [Abstract] [CrossRef] [Google Scholar]

Citations & impact 


Impact metrics

Jump to Citations

Citations of article over time

Alternative metrics

Altmetric item for https://www.altmetric.com/details/146327318
Altmetric
Discover the attention surrounding your research
https://www.altmetric.com/details/146327318

Article citations

Similar Articles 


To arrive at the top five similar articles we use a word-weighted algorithm to compare words from the Title and Abstract of each citation.


Funding 


Funders who supported this work.

NICHD NIH HHS (1)