Modeling of the Hemodynamic Responses in Block Design fMRI Studies

Simulation of Functional Magnetic Resonance Imaging Signal Time Course

Functional magnetic resonance imaging signal time courses were simulated using known HRFs to test the performance of different HRF models. Three ground truth HRFs were created (Supplementary Figure S1). In Simulation 1, the HRF was generated using the balloon model.^{17, 18} The HRF was generated using a single input at 1.5 seconds with 0.1-second duration and using typical parameters obtained empirically in a previous study: neuronal efficacy 0.54, signal decay 1.54, autoregulation 2.46, transit time 0.98, stiffness parameter 0.33, resting oxygen extraction 0.34, and resting blood volume fraction 0.06. ¹⁸ The HRF in Simulation 2 was generated using the sum of four inverse logit functions to create four segments of HRF time course:

h ( t ) = ∑ i = 1 4 A i 1 1 + e ( t − T i ) ∕ D i ,

(6)

in which A ₁ =−0.2, A ₂ =1.8, A ₃ =−1.8, A ₄ =0.2, T ₁ =0.1, T ₂ =4, T ₃ =10, T ₄ =20, D ₁ =0.8, D ₂ =1, D ₂ =1, and D ₃ =1.2. These parameters were determined empirically as generating an HRF shape with summary features (H, W, T, and O) in keeping with published values.^{19, 20} In Simulation 3, we used the sum of three gamma functions to generate the HRF:

h ( t ) = ∑ i = 1 3 ( A i t α i − 1 β i α i e − β i t Γ ( α i ) ) ,

(7)

in which A ₁ =−0.2, A ₂ =10, A ₃ =−3.6, α₁=1.5, α₂=6.6, α₃=15, β₁= 0.8, β₂=0.8, and β₃=1. Figure 1 illustrates the procedures we used to simulate the fMRI time series. We used a block design stimulus function comprising 16 alternating blocks of rest and active conditions (8 for each condition). Each active block included 16 trials of 0.2-second duration with 0.8 seconds gap (Figure 1B). Random noise was added to achieve SNR levels of 50, 100, and 150 (Figure 1D). Noise was added without consideration of autocorrelation:

T C = T C × ( 1 + σ × r a n d n ) ,

(8)

where TC is the time course of simulated fMRI data; σ is the s.d. (0.02, 0.01, and 0.067 for SNR of 50, 100, and 150, respectively); and the randn refers to a vector of pseudorandom values drawn from the standard normal distribution with the same length as TC. The range of SNR levels was selected given previous studies showing that the minimum SNR for detecting signal change in fMRI studies is 69 and the highest possible SNR is 154. ²¹ The simulated signal was sampled according to a typical fMRI acquisition with TR=2.1 seconds (Figure 1E). Each simulation was run 100 times.

Hemodynamic Response Function Model Selection

Hemodynamic response function modeling was performed using an in-house MATLAB (MATLAB R2010a, The MathWorks, Natick, MA, USA) toolkit, sHRF, which is publicly available at our website (http://www.cai.uq.edu.au/shrf-toolkit). The fMRI time courses were high-pass filtered with a 128 seconds discrete cosine basis set and then modeled as the convolution of the HRF (Table 1) and the stimulus function: ²²

T C = h ( t ) ∗ u ( t ) + C ,

(8)

where h(t) and u(t) are the HRF and stimulus function, respectively. For the fitting we used a modified constrained Nelder–Mead Simplex algorithm that allows the use of constraints specified as parameter bounds (http://www.mathworks.com/matlabcentral/fileexchange/8277). Fitting was initialized using parameters, which generated an initial HRF shape that closely matched that of the canonical HRF. Initial values and constraints for each model are summarized in Table 1. Parameter constraints were set empirically to cover plausible HRFs with the following characteristics: an initial dip with a minimum at 0 to 3 seconds of magnitude 0% to 5%; peak BOLD signal at 2 to 10 seconds of magnitude 0% to 15%; signal undershoot with a minimum at 6 to 25 seconds of magnitude 0% to 10%. The five models cover the same space of possible HRFs except where the functional form imposes a limitation in the ability to describe the HRF. For example, Model I covers the same space for the amplitude of the BOLD peak as other models but does not allow temporal variation; Model II covers the same space for the BOLD peak and undershoot as Models III to V but the functional form does not allow an initial dip. The parameter C allows baseline adjustment of the fitted time course to match that of the simulated time series (Figure 1D). The trueness of estimated HRF parameters or the closeness of agreement between measured and ground truth values was assessed by calculating the relative deviation from the values for the HRF used to generate the simulated fMRI time course, i.e. the difference between estimated and ground truth values divided by the ground truth value. The measure of precision, or relative uncertainty, used in this study was the interquartile range of relative deviation from the values for the HRF used to generate the simulated fMRI time course. Trueness and precision were summarized using boxplots of the simulation results created in MATLAB (MATLAB R2010a, The MathWorks). The identifiability of HRF models was assessed using their Akaike weights, ²³ which is the probability of each model transformed from the Akaike Information Criterion (AIC):

w i = e − 1 2 Δ i ( A I C ) ∕ ∑ j = 1 M e − 1 2 Δ j ( A I C ) , Δ i ( A I C ) = A I C i − min ( A I C ) ,

(9)

in which M is number of models tested and Δ _i (AIC) is the difference between the AIC of each model and that of model with the smallest AIC value. In this study, we used the AIC for finite samples: ²⁴

A I C = n ln ⁡ ( R S S ∕ n ) + 2 k ( k + 1 ) ∕ ( n − k − 1 ) ,

(10)

where RSS is the residual sum of squares between the fitted and fMRI time courses, n is the sample size, and k is the number of parameters. The HRF model with the highest AIC weight was selected for subsequent analyses of parameter estimation and reproducibility.

Parameter Estimation from Block Design Functional Magnetic Resonance Imaging Data

The noise in measured fMRI time series and sparse sampling in acquired data may lead to HRF parameter mis-specification. To examine this, we used the balloon model simulation^{17, 18} as the ground truth. The HRF parameters were varied in three ways: (1) H only, (2) T and O only, and (3) W, T, and O only. As illustrated in Supplementary Figure S2, varying T unavoidably changed O and varying W unavoidably changed T and O. The temporal SNR was set at 100, a typical SNR for fMRI data according to previous reports.^{25, 26} The highest variability across subjects observed in the literature is 25%.^{1, 3, 19} Hence, we varied ground truth HRF parameter values by 25%, 15%, and 10%. The simulation procedure described above in the section on HRF model selection was employed, although here only Model III was used to estimate the HRF. The simulation was executed 100 times for each type and level of parameter variation. The estimated HRF parameters were compared using the paired t-test in SPSS20 (IBM, Madison Avenue, New York City, NY, USA).

Variations in the Block Design Paradigm

To investigate the influence of variations in the block design paradigm on the trueness and precision of HRF modeling, the experimental conditions of the block design were systematically varied: (1) the number of blocks was varied from 1 to 8 with one block increments; (2) block length was varied from 4 to 16 seconds with 2-second increments; (3) stimulus duration was varied from 0.2 to 6.2 seconds with 1-second increments and a stimulus duration of 16 seconds in individual blocks (long block design) was also evaluated; (4) the gap between individual stimuli in each block was varied from 0.8 to 6.8 seconds with 1-seconds increments. Simulation of the fMRI time course was carried out as described above using the balloon model to generate ground truth, a SNR of 100, and Model III to estimate the HRF. Each simulation was executed 100 times. Trueness and precision for each paradigm were calculated as described above.

Experimental Functional Magnetic Resonance Imaging Time Course Data

The data we analyzed were acquired as part of a prior fMRI study, ¹¹ the Queensland Twin Imaging Study. ²⁷ Retrospective utilization of the data was approved by the Research Ethics Committee of the Queensland Institute of Medical Research and The University of Queensland in compliance with the Australian National Statement on Ethical Conduct in Human Research. Functional MRI data from 30 male and 52 female healthy young adult twins of mean age 22.5±2.5 (s.d.) years were included. For all participants, the fMRI scans were repeated 2 to 6 months later (mean 117±56 days).

Participants performed the 0- and 2-back versions of the N-back working memory task. The detailed fMRI experimental procedure is described in a previous report. ²⁸ In the N-back task, a series of numbers is presented on a screen. The 0-back condition required participants to respond to the number currently shown on the screen. The 2-back condition required participants to respond to the number presented 2 trials earlier. The number was presented for 200 milliseconds with an 800-millisecond interval between stimuli and 16 trials per block. In total, 16 alternating blocks were performed for the two conditions (eight blocks per condition).

The three-dimensional T1-weighted MR image and echo planar imaging data were acquired on a 4T Bruker Medspec whole body scanner (Bruker, Hanau, Germany). Three-dimensionalT1-weighted images were acquired using an MP-RAGE pulse sequence (TR=2500 milliseconds, TE=3.83 milliseconds, T1=1500 milliseconds, flip angle=15 degrees, 0.89 × 0.89 × 0.89 mm). For each participant, 127 sets of echo planar imaging data (TR=2.1 seconds, TE=30 milliseconds, flip angle=90, 3.6 × 3.6 × 3.0 mm) were acquired continuously during the tasks.

The fMRI data were analyzed using Statistical Parametric Mapping (SPM8, the Wellcome Trust Centre for Neuroimaging, London, UK). The first five echo planar imaging volumes were discarded to ensure that tissue magnetization had reached steady state. The spatial preprocessing included 2-pass motion correction ²⁹ and spatial normalization to the average brain T1 template ³⁰ implemented in SPM8. Normalized volumes were smoothed with an 8 × 8 × 8 mm full width half maximum Gaussian kernel. In a first-level analysis of individual subject data, we determined locations of activation using the general linear model with a finite impulse response basis function. The 2-back minus 0-back contrast images were then entered into a group level, random-effects one-sample t-test to identify the common activation voxels (P<0.05 with family wise error rate adjustment for multiple comparisons).

Four regions identified at the group level were selected for HRF modeling: left and right middle frontal gyrus and left and right angular gyrus (Supplementary Figure S3). Volumes of interest were defined as the overlap of an existing probabilistic atlas of each structure in stereotaxic coordinate space ³¹ with group activation regions. For each participant, the fMRI time course was extracted by averaging the signal intensity at each time point in the voxels with the top 12.5% of SPM t statistics within each VOI. The top 12.5% is used here to accommodate individual variation in the pattern of functional activation within each VOI. The threshold was selected empirically as the value when the extracted time course stabilized as the percentage was decreased from 100, 50, 25, to 12.5. The HRF model selected from the simulation study was fitted to the extracted time series and used to estimate HRF parameters. Akaike Information Criterion weights were also calculated for each HRF model using real time series data.

Analysis of Test–Retest Reliability

For each brain region and each HRF parameter, the intraclass correlations of the two experimental sessions were calculated using SPSS20 (IBM) using a two-factor mixed effects model and tests of significant difference in ICC from zero performed. ³² Because the participants were biologically related monozygotic and dizygotic twins, the intraclass correlations of the HRF heights from the first-born participant in each twin pair were also evaluated to assess whether relatedness between participants affected the reliability results. The reliability of the performance measure used in the experiment has been reported previously. ¹¹

RESULTS

Simulation Study

Figure 2 summarizes the HRF parameter estimates from the simulations. Parameter estimation was considered accurate if estimated values did not differ significantly (i.e. P>0.05) from the corresponding ground truth values. Height was estimated accurately for all SNR levels and all models except for Model V at an SNR of 50 in simulation 1 (P<0.05, one sample t-test). For Model I, T, W, and O are fixed and there is only one free parameter. Models II, III, and IV were able to estimate T, W, and O of the HRF for all noise levels accurately. Model V was able to estimate T and W accurately at all noise levels but did not estimate O accurately (P<0.05 for all three simulations). The central line and lower and upper box boundaries of the boxplots in Figure 3 respectively represent the median, 25th and 75th percentiles of 100 simulations. Models with more free parameters showed a lower precision (wider range) of parameter estimation but higher trueness (smaller relative deviation) than those with fewer parameters. Model I had the least variation but the highest bias of parameter estimation. For all three simulations, the variation and bias of parameter estimation decreased with increasing SNR for all models. For Models II to V, T was estimated with the highest precision and trueness followed by W, H, and O.

OPEN IN VIEWER

Figure 2.

Box-and-whisker plot of the relative differences of estimated hemodynamic response function (HRF) parameters from the true values in the simulation study. Outliers were defined as exceeding three times the s.d.; all data including outliers were included in the analysis. Rows from top to bottom are the results for height (H), time to peak (T), width (W), and onset (O). Columns from left to right are results for Simulation 1 (‘balloon’ model), Simulation 2 (sum of four inverse logit functions), and Simulation 3 (sum of three gamma functions).

OPEN IN VIEWER

Figure 3.

Boxplots of estimated hemodynamic response function (HRF) features with systematic variation in attributes of the block design paradigm. The columns from left to right relate to estimated height, time to peak, width, and onset, respectively. The rows from top to bottom relate to variation in the number of blocks (A), the length of blocks (B), the duration of stimuli in each block (C), and the gap between individual stimuli (D). Outliers were defined as exceeding three times the s.d.. Outliers were not excluded from data analysis.

Table 2 summarizes the AIC of different models for data with different noise levels. Model III had the highest AIC weights and Model I the lowest across different simulations and noise levels.

OPEN IN VIEWER

Table 2.

Akaike weights (W _i ) for HRF models for computer simulations and real fMRI data ^a

HRF model	Simulation 1			Simulation 2			Simulation 3			Real data
	150	100	50	150	100	50	150	100	50
III (2 Gamma functions 6 parameters)	0.34	0.31	0.25	0.36	0.33	0.27	0.4	0.33	0.3	0.41
IV (3 Gamma functions 9 parameters)	0.27	0.23	0.2	0.34	0.22	0.22	0.32	0.25	0.22	0.36
II (2 Gamma functions 5 parameters)	0.22	0.21	0.23	0.25	0.3	0.23	0.26	0.26	0.23	0.16
V (3 Inverse logit functions 7 parameters)	0.16	0.2	0.2	0.05	0.15	0.21	0.02	0.15	0.18	0.07
I (Canonical function)	0.01	0.05	0.12	<0.01	<0.01	0.07	<0.01	0.01	0.07	<0.01

fMRI, functional magnetic resonance imaging; HRF, hemodynamic response function; SNR, signal-to-noise ratio.

a

The Akaike weights were calculated for each model at a specific SNR level (150, 100, and 50) with different ground truth generation method (Simulation 1, 2, and 3) on simulation data and on real fMRI data. The simulation was run 100 times for each condition.

Results for HRF parameter estimation are summarized in Table 3. The selected HRF model (Model III) was able to identify parameter differences exceeding 10% with SNR of 100. There were no instances in which significant changes in one parameter resulted from variation in another parameter.

OPEN IN VIEWER

Table 3.

Estimation of HRF parameter variation ^a

		Estimated parameters (±s.d.) with variations in percentage
P	Δ=0	ΔH=10	ΔH=15	ΔH=25	ΔT=10O=18	ΔT=15O=32	ΔT=25O=50	ΔW=10ΔT=9	ΔW=15ΔT=13O=4	ΔW=25ΔT=20O=4
H (%)	1.12±0.46	1.21±0.49	1.31± 0.32∗	1.41±0.38∗∗	1.14±0.39	1.06± 0.39	1.11±0.37	1.13± 0.45	1.17±0.39	1.18±0.35
T (s)	4.70±0.92	4.65±0.91	4.71±0.78	4.45±0.65	5.08±0.81∗	5.53±0.75∗∗	5.82± 0.68∗∗	5.14±0.93∗	5.14±0.84∗	5.35±0.85∗∗
W (s)	3.85±1.07	4.05±1.07	4.10±0.96	4.00±0.98	4.08±1.01	4.23±1.25	4.15±0.73	4.53±1.09∗∗	4.51±1.01∗∗	4.77±0.97∗
O (s)	1.90±0.88	1.79±0.84	1.83±0.74	1.74±0.82	2.20±0.90	2.41± 0.85∗	2.79±0.89∗	1.98±0.81	2.02±0.89	2.08±0.89∗

HRF, hemodynamic response function; H, height; O, onset; P, HRF parameters; T, time to peak; W, width; Δ, variation (no variation is introduced if Δ=0 or the parameter is absent).

a

The ∗ and ∗∗ denote P<0.05 and P<0.001, respectively.

Results of parameter estimation with different block design paradigms are summarized in Figure 3 using boxplots constructed in the same manner as Figure 2. For less than six blocks, mean estimated H, W, and O generally had relative errors greater than 0.1 and an interquartile range greater than 0.5 (Figure 3A). Trueness was similar across paradigms with different block lengths. However, precision was increased with reduced block length (Figure 3B). Both trueness and precision of parameter estimation increased with greater stimulus duration within each block (Figure 3C). Trueness and precision also increased when the gap between stimuli was increased to 3.8 seconds and remained similar with larger gaps (Figure 3D).

Real Functional Magnetic Resonance Imaging Data

Akaike Information Criterion weights for the models using real fMRI data were consistent with results from simulations (Table 2). Model III had the highest AIC weights. Activation maps between the two experimental sessions were highly reproducible (Supplementary Figure S3). Quantitative assessment of the repeatability of the activation map between two experimental sessions has been reported previously. ¹¹ The ICC(3,1) of HRF parameters estimated in the two experimental sessions using Model III are summarized in Table 4. Height and T are highly reproducible, W is moderately reproducible, and reproducibility of O is low. The 95% confidence intervals of all parameters overlap between the two experimental sessions. Height is lower in session 2 than in session 1. A practice effect was also observed in the performance data. ¹¹ Biologic relatedness between participants did not affect reliability results. The intraclass correlations of HRF heights from the first-born participant in each twin pair (N=41) for left and right middle frontal gyrus and left and right angular gyri were 0.64 (P<0.01), 0.58 (P<0.01), 0.63 (P<0.01), and 0.64 (P<0.01), respectively.

OPEN IN VIEWER

Table 4.

The ICC(3,1) of HRF parameters between the two experimental sessions ^a

Structures	Parameters
	Height	Time to peak	Width	Onsets
L-MFG	0.6∗∗ (1.23, 1.56; 1.03, 1.26)	0.53∗∗ (3.57, 4.09; 3.73, 4.21)	0.56∗∗ (3.49, 4.08; 3.84, 4.27)	0.22 (1.17, 1.5; 1.21, 1.49)
R-MFG	0.52∗∗ (1.21, 1.52; 1.05, 1.28)	0.5∗ (3.7, 4.29; 3.94, 4.52)	0.45∗ (3.75, 4.39; 4.26, 5.08)	0 (1.22, 1.54; 1.23, 1.59)
L-AG	0.74∗∗ (1.28, 1.56; 1.02, 1.29)	0.7∗∗ (3.66, 4.19; 3.78, 4.35)	0.48∗ (3.75, 4.27; 3.91, 4.74)	0.25 (1.18. 1.49; 1.24, 1.57)
R-AG	0.6∗∗ (1.12, 1.37; 1.02, 1.26)	0.45∗∗ (3.68, 4.29; 3.8, 4.43)	0.33∗ (4.02, 4.74; 4.28, 4.91)	0.42∗ (1.12, 1.47; 1.15, 1.59)

AG, angular gyrus; HRF, hemodynamic response function; L, left, MFG: middle frontal gyrus, R, right.

a

∗P0.05 and ∗∗P<0.001. The numbers in the parentheses are the 95% confidence interval of the parameters (session 1 lower bound, session 1 upper bound; session 2 lower bound, session 2 upper bound).

DISCUSSION

Modeling of the HRF using fMRI is a noninvasive method to investigate the brain's hemodynamic response to neuronal activation. In this study, we found that with block design fMRI paradigms, most HRF models were able to estimate HRF parameters except for Model V that had difficulty estimating onsets. The sum of two gamma functions with six free parameters had the greatest identifiability of the five models tested but the accuracy of parameter estimation varied with block design attributes. H, T, and W were reproducible between two experimental sessions.

We used computer simulations with a known ground truth to assess the performance of the HRF models. In the simulations, Gaussian noise was added at each SNR level without consideration of autocorrelation. However, fMRI noise typically exhibits temporal dependence (autocorrelation); therefore, the SNR level is slightly lower than in practice. Parrish et al ²¹ demonstrated that a minimum SNR of 69 is required to detect a 1% BOLD signal change with 112 image volumes. Our results were similar, showing that Models I to IV could estimate H accurately with a SNR of 50 and 122 image volumes. The simulation results for the five models showed that models with more free parameters had lower precision but higher trueness of HRF parameter estimation (Figure 3). The canonical HRF (Model I) had the least variance, but the estimated height was more likely to be biased. The tradeoff between the number of parameters and precision and trueness motivated the use of AIC weights for model selection.

A key consideration in HRF model selection is the number and nature of a priori assumptions. Recently, Lindquist et al ⁴ assessed the performance of different HRF models for event-related designs in terms of bias, mis-specification, and statistical power. Seven models were compared: the canonical HRF, canonical HRF with temporal derivatives, canonical HRF with temporal and dispersion derivatives, finite impulse response, smoothed finite impulse response, sum of two gamma functions, and the sum of three inverse logit functions. The study suggested the sum of two gamma functions is not optimal for fitting noisy data from event-related studies and should only be used on regions where it is known that there is signal present. ⁴ Our study focused on HRF modeling using block design fMRI data in which signal is present. As illustrated in the simulation study (Figure 3), models with few assumptions are more flexible or may handle HRF shapes with unexpected response behaviors more accurately than models with many fixed parameters. However, a greater number of free parameters may also lead to overfitting of noise, lowering precision. This may explain the relatively poor performance of Model V in our study. In Lindquist et al, ⁴ superior performance for Model V compared with gamma function models for event-related designs that have a high HRF estimation efficiency ⁸ and high temporal resolution (TR=0.5 seconds),^{4, 15} most likely reflects Model V's greater flexibility.^{4, 15} In contrast, in the present study of block designs with low estimation efficiency, prior knowledge implicitly encoded in the gamma function models constrained them from overfitting noisy data. We used AIC weights to judge model identifiability. The AIC is a measure of the relative goodness-of-fit of a statistical model compared with others, with a penalty given for using extra parameters. ³³ The AIC with a correction for finite sample sizes was used in this study (equation (10)); the correction is minimized when the sample size is large. ²⁴ Based on AIC values, Wagenmakers et al ²³ further proposed Akaike weights, which may be directly interpreted as conditional probabilities for each model. For example, an AIC weight of 0.41 for Model III with real data in our study indicates that this model has a 41% probability of being the best model of those tested. The AIC values of different models are consistent between the simulations and real fMRI data. The results suggest that the sum of two gamma functions with six parameters is the best HRF model when it is known that there is signal present. With a typical fMRI SNR of 100, a two-gamma-function, six-parameter model can discern 15% HRF feature variation without misattributing variation in one parameter to another. Inter-and intra-subject variability of the HRF has been studied previously in event-related fMRI experiments.^{1, 3, 19, 34, 35, 36} In the present examination of a block design study, H, T, and W were reproducible between two experimental sessions, supporting the validity of these HRF parameters as biologic measures. The repeatability results may also provide insight into neural responses and neurovascular coupling. For example, the reproducibility of H between trials implies that the neural firing rate as well as the vascular response induced by the same task across different sessions for each individual was very similar. In previous studies,^{5, 37, 38} the height of the BOLD response to visual stimulation was negatively correlated with resting GABA concentration in the visual cortex. If this relationship is widespread in the brain, our results suggest that endogenous factors, such as neurotransmitter concentration, which regulate the balance between excitation and inhibition in the individual brain may also remain relatively stable. The reproducibility of T and low reproducibility of O has implications for the ability to infer synchrony and directional influence among neural populations during task performance and in the resting state from block design fMRI data. Our findings on HRF test–retest reliability, particularly for temporal parameters, are in broad agreement with results from studies optimized for HRF estimation by the use of event-related designs and long baseline periods,^{1, 36} in which time to peak had the greatest test–retest reliability and time when the BOLD response begins to rise steeply from baseline (similar to O in this study) had the largest test–retest variation.

In this study, we represented neural activity by the stimulus time course, an assumption underlying commonly used methods for statistical parametric mapping. To the extent that this assumption is violated, it is likely that the HRF model parameters estimated here contain a component of neural activity in addition to vascular reactivity. Several techniques have been recently proposed to dynamically filter fMRI time courses to estimate neural activity as well as the vascular response explicitly.^{39, 40} Application of these techniques in future studies may allow the reproducibility of the HRF to be estimated more specifically.

Many HRF models have been proposed in the literature. Our study was restricted to five models that are commonly used or have been shown in prior studies to be less susceptible to mis-specification errors (sum of inverse logit functions) or parameter bias (sum of two gamma functions). In this study, we examined vascular responses to neural activity evoked by a single cognitive function in four brain structures that demonstrated robust activation. Further studies using other activation tasks are needed to assess generalizability of our findings to other brain regions and brain functions.

Block designs are not optimal for HRF estimation and random event-related designs and m-sequence designs have been shown to have higher estimation efficiency. ⁹ Nonetheless, our data on model selection and on the reliability of HRF estimation should increase the utility of existing block design fMRI data. Indeed, the methodology presented in this paper will be used to analyze the heritability of HRF parameters in a large twin cohort using block design fMRI data that were originally collected to address heritability of spatial patterns of functional activation with a working memory task. ¹¹ In our tests of the effects of varying the attributes of the block design paradigm, we found that at least six blocks are required to model the HRF accurately. The simulations showed that prolonging the length of blocks does not affect the trueness of the parameter estimation but decreases precision. This result is in keeping with HRF modeling being more sensitive to the initial increase and the final decrease of the signal than the plateau stages in the middle of block. Prolonging the block introduces additional noise reducing precision. Decreasing the frequency of stimuli in each block by increasing the gap between stimuli increases the accuracy of parameter estimation. This finding agrees with a previous theoretical model describing the tradeoff between estimation efficiency and detection power ⁸ in the sense that block designs with decreasing stimulus frequency in each block bear a greater resemblance to event-related designs.

In conclusion, we investigated trueness, precision, model identifiability, and test–retest reproducibility of parameter estimates for HRF models with a block design paradigm using simulated and real fMRI data. The HRF features of height, time to peak, and width were reproducible between test sessions and may be useful as measures to characterize the coupled vascular response to neural activity in individual subjects.