Fast fMRI responses supported by inter-voxel variability of hemodynamic response functions
Jingyuan E. Chen1,2, Jonathan R. Polimeni1,2,3, Nina E. Fultz1, Gary Glover4, and Laura D. Lewis1,5

1Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Boston, MA, United States, 2Department of Radiology, Harvard Medical School, Boston, MA, United States, 3Harvard-Massachusetts Institute of Technology Division of Health Sciences and Technology, Cambridge, MA, United States, 4Radiology, Stanford University, Palo Alto, CA, United States, 5Biomedical Engineering, Boston University, Boston, MA, United States


In this abstract, we use very high-resolution fMRI to evaluate inter-voxel variability of hemodynamic response functions (HRF) elicited by a visual task. We show that HRFs exhibit wide variability across voxels, with a subset of HRFs occurring at temporal scales much faster than the canonical HRF, which we consider as one mechanism that underlie fast neural-activity-related dynamics. Influence of spatial resolution on HRF dynamics are also evaluated, we show that the use of small voxels (without spatial averaging) sampling parenchymal signals may enable faster observed hemodynamics.


Recent resting state and task-driven fMRI studies have demonstrated detectable neural fluctuations at high frequencies (>0.2 Hz) that are far larger than predicted by the canonical hemodynamic response function (HRF)[1,2]. Without contradicting the classical hemodynamic response models, several possible mechanisms have been proposed, attributing fast fMRI dynamics to well-known nonlinear hemodynamic changes elicited by short-duration stimuli[3], altered hemodynamics during low-amplitude neural activation[4], or as of yet un-identified non-BOLD mechanisms[4]. Here, we show using very high resolution fMRI that as HRFs exhibit wide variability even within small cortical regions, a considerable number of individual voxels demonstrate much faster responses than the grand average, which resembles the canonical HRF. Thus, the use of small voxels (without spatial averaging) sampling parenchymal signals may enable faster observed hemodynamics .


Acquisition & Experiments: five healthy volunteers participated in this experiment after providing written informed consent, with four participants scanned at 3T (Siemens Prisma scanner, TR/TE/FA = 780ms/34ms/53o, voxel size = 1.23×1.23×1.5 mm3, 9 coronal slices covering the calcarine sulcus were acquired) and one subject scanned at 7T (a whole-body scanner (Siemens Healthineers, Erlangen, Germany), TR/TE/FA = 1010ms/28ms/70o, voxel size = 0.8 mm iso. in-plane acceleration factor = 5,15 coronal slices covering calcarine sulcus were acquired). Each subject underwent 5~6 event-related visual task scans (12 Hz flickering checkerboard pattern) with stimulus presentation following a pseudo-random m-sequence paradigm (with base-interval being 0.5 s, and 256 s-long total duration per scan). Voxel-wise HRFs: after motion/slice timing correction and inter-scan registration, we resampled each voxel’s time course to match the paradigm of m-sequence (sampling rate = 0.5 s), then de-convolved the impulse responses triggered by the 0.5 s long visual trial for each voxel. Smooth HRF bases were generated in MATLAB following the FSL FLOBS concept[5], with the parameters outlined in Fig. 1A. Two optimal basis sets comprising 3 and 6 basis elements (Fig. 1B) were employed to fit voxel-wise task impulse responses, termed as ‘HRF_flob3’ and ‘HRF_flob6’ respectively. Inter-voxel variability of HRF temporal characteristics (timing and frequency responses): metrics including time-to-peak (TTP), full width at half maximum (FWHM), and normalized frequency responses (normalized by the frequency amplitude at 0.1 Hz to enable comparisons across voxels), were computed for voxels with T-scores > 4. Influence of voxel size on HRF timing: the 0.8 mm isotropic resolution 7T data were spatially smoothed with varying sizes of Gaussian kernels (FWHM = 1/2/4 mm) to evaluate how spatial pooling influences the timing and temporal frequency content of observed HRFs.

Results & Discussion

Robust activation surrounding the calcarine sulcus was observed in all five subjects—results from the 7T subject and one representative 3T subject are shown in Fig. 2. As shown in both Fig. 3 and 4 (‘no smoothing’), remarkably variable TTPs and FWHMs across task-active voxels were observed for each subject, with the slowest TTPs/FWHMs being approximately twice long as the fastest/narrowest ones. As a result, normalized frequency responses at higher frequencies also exhibit prominent variability—with the maximum response being ~10-fold larger than the minimal response at 0.4 Hz. A large number of voxels exhibit much faster responses than the canonical HRF[6] (TTP = 5 s, FWHM = 5.3 s, normalized frequency amplitude of 0.17 at 0.2 Hz and 0.008 at 0.4 Hz).

As expected, coarser effective voxel sizes imposed by spatial smoothing reduces the dispersion of HRF timings (Fig. 4A, compare across rows), and obscures the fastest hemodynamic responses. Notably, the mean impulse responses appear faster following smoothing (cf. the center distribution across rows in Fig. 4A, Fig. 4B). This potentially counter-intuitive finding can be explained by considering that the slowest HRFs are expected to be found in large veins, which are sparse in the cortex, therefore averaging parenchymal signals with large vein signals will result in a net increase in HRF dynamics.


Here, by analyzing HRF variability at fine spatial scales, we show that HRFs may give rise to frequency responses much faster than those predicted by the canonical model. Not only does this suggest that the observed dynamics will be a function of spatial resolution, but it also suggests a strategy for increasing the speed of the observed responses by using small voxels confined to the parenchyma, while avoiding large draining vessels on the surface which exhibit delayed, downstream responses. Conversely, low spatial resolution reduces HRF variability and limiting the observable frequency upper-bound, while causing an increase in HRF speed in regions exhibiting significant percent signal change.


This work was supported in part by the NIH NIBIB (grants P41-EB015896, and R01-EB019437), NINDS (R21-NS106706) by the BRAIN Initiative (NIH NIMH grants R01-MH111438 and R01-MH111419), and by the MGH/HST Athinoula A. Martinos Center for Biomedical Imaging; and was made possible by the resources provided by NIH Shared Instrumentation Grants S10-RR023043 and S10-RR019371.


[1] Lee, H.L., Zahneisen, B., Hugger, T., Levan, P., Hennig, J., 2013. Tracking dynamic resting-state networks at higher frequencies using MR-encephalography. Neuroimage 65, 216-222.

[2] Lewis, L.D., Setsompop, K., Rosen, B.R., Polimeni, J.R., 2016. Fast fMRI can detect oscillatory neural activity in humans. Proc Natl Acad Sci U S A 113, E6679-E6685.

[3] Lewis, L.D., Setsompop, K., Rosen, B.R., Polimeni, J.R., 2018. Stimulus-dependent hemodynamic response timing across the human subcortical-cortical visual pathway identified through high spatiotemporal resolution 7T fMRI. Neuroimage 181, 279-291.

[4] Chen, J.E., Glover, G.H., 2015. BOLD fractional contribution to resting-state functional connectivity above 0.1 Hz. Neuroimage 107, 207-218.

[5] Woolrich, M.W., Behrens, T.E.J., Smith, S.M. 2004. Constrained linear basis sets for HRF modelling using Variational Bayes. NeuroImage 21, 1748-1761.

[6] https://www.fil.ion.ucl.ac.uk/spm/software/spm12/.


Fig. 1 A: Parameters defining the space spanned by the HRF basis. (B) Sets of individual temporal basis elements chosen to estimate the HRFs within individual voxels.

Fig. 2: T-score map of visual task activation (averaged across all task scans) in the 7T subject and one representative 3T subject. The statistical maps are overlaid on top of raw BOLD-weighted images (gradient-echo EPI).

Fig. 3 Histograms of HRF timing and frequency responses across voxels (with T-score > 4) in four 3T subjects, without spatial smoothing. The distribution of HRF waveform parameters is large, causing a frequency response distribution spanning nearly an order of magnitude.

Fig. 4 Influence of spatial resolution (smoothing) on HRF dynamics, illustrated using the 0.8 mm iso. res. 7T dataset. A: Plotting the distribution of HRF timing and frequency responses across voxels shows that smoothing reduces the range of observed HRF waveforms (the identical set of voxels were compared across different smoothing kernels). B: Illustration using single-voxel time-series data that lower spatial resolution may counter-intuitively result in faster hemodynamics (seen in a shifting of the distribution towards shorter TTP and FWHM values) in those regions demonstrating significant task effects.

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)