Towards whole brain mapping of the haemodynamic response function

Abstract

Functional magnetic resonance imaging time-series are conventionally processed by linear modelling the evoked response as the convolution of the experimental conditions with a stereotyped haemodynamic response function (HRF). However, the neural signal in response to a stimulus can vary according to task, brain region, and subject-specific conditions. Moreover, HRF shape has been suggested to carry physiological information. The BOLD signal across a range of sensorial and cognitive tasks was fitted using a sine series expansion, and modelled signals were deconvolved, thus giving rise to a task-specific deconvolved HRF (dHRF), which was characterized in terms of amplitude, latency, time-to-peak and full-width at half maximum for each task. We found that the BOLD response shape changes not only across activated regions and tasks, but also across subjects despite the age homogeneity of the cohort. Largest variabilities were observed in mean amplitude and latency across tasks and regions, while time-to-peak and full width at half maximum were relatively more consistent. Additionally, the dHRF was found to deviate from canonicity in several brain regions. Our results suggest that the choice of a standard, uniform HRF may be not optimal for all fMRI analyses and may lead to model misspecifications and statistical bias.

Keywords: BOLD response; HCP; HRF; fMRI; haemodynamic.

Figure 1.

Average fitted responses and dHRFs across subjects. Task-specific mean fitted BOLD responses and mean estimated dHRFs for M1 (top, a and b) and M2 (bottom, c and d). Data are averaged in each responding area and then across subjects. In all plots, the average response is represented with a black solid line while the coloured shades around the mean represent the 95% confidence interval of the mean across subjects.

Figure 2.

Across-subject variability in amplitude and timing parameters. Across-subjects distribution of amplitude (a) and timing parameters (latency, TTP, and FWHM) for M1 (top, a and b) and M2 (bottom, c and d). Parameters were extracted from the relevant model and calculated for each task and subject using the spatially averaged task-based response and dHRF. The box plots show the mean (open circle), the median (horizontal line), and the interquartile range (box extremes); the whiskers extend to the data point closest to 1.5 times the interquartile range, above or below the first and third quartiles, respectively; grey dots represents outlier points ranging beyond these thresholds.

Figure 3.

Correlations between M1 parameters. (a) Correlation between TTP and latency (significant, R = 0.16, p = 0.00024). (b) Correlation between TTP and FWHM (significant, R = 0.35, p < 0.0001) and (c) correlation between latency and FWHM (not significant, R = −0.08, p = 0.095).

Figure 4.

Adjusted R² of M2 vs M1. Plot of adjusted R² for all subjects and tasks. Adjusted R² of regression model M1 was always significantly superior to that of model M2 (p < 10⁻⁸ for all tasks).

Figure 5.

Canonicity of the dHRF. Voxelwise similarity, expressed as z-score, of the dHRF to the SPM12 canonical HRF (a), the relevant SD (b) and sagittal section and scale (c). For each voxel and subject, the task-averaged dHRF was obtained by averaging the dHRF shapes of all tasks for which the voxel was responding and was then correlated to the canonical HRF. Regions showing the maximal canonicity are in the occipital cortex. Several regions, mainly in the insular somatosensory and motor cortex and frontal pole showed some asymmetry between hemispheres.

Figure 6.

Whole brain maps of the investigated BOLD parameters, estimated with M1 model. In each row, from left: axial map of the parameter, associated standard deviation across subjects, sagittal section and scale. From top to bottom: (a)–(c) amplitude of the fitted Continued.BOLD response, expressed as percentage variation of BOLD signal over the baseline. The amplitude of the fitted response varied appreciably across the cortex. In particular, a larger BOLD signal change was observed in sensory areas, such as the primary and supplementary motor cortex, intracalcarine and supracalcarine cortex, and in the caudate nucleus and cingulate gyrus. (d)–(f) mean latency of the fitted BOLD response. Latency appeared moderate and reproducible across subjects in occipital regions, while some temporal areas were characterized by slower responses. (g)–(i) Time To Peak of the deconvolved HRF. TTP showed a high spatial variability at the voxel level. Higher TTP values were found in the posterior division of the temporal gyrus while lower values were found in the supplementary motor cortex. (j)–(l) Full Width at Half Maximum of the dHRF. FWHM was found to be higher in the frontal medial cortex, lateral occipital cortex and angular gyrus, motor cortex (precentral gyrus), cingulate gyrus, and postcentral gyrus. For each parameter, the voxelwise magnitude was first averaged across the responding tasks and then the mean and SD across subjects were calculated.