NEUROIMAGE 8, 360–369 (1998)

ARTICLE NO. NI980369

The Variability of Human, BOLD Hemodynamic Responses

G. K. Aguirre, E. Zarahn, and M. D’Esposito
Department of Neurology, Hospital of the University of Pennsylvania, Philadelphia, Pennsylvania 19104-4283

Received November 3, 1997

that follows a brief period of neural activity can be

Cerebral hemodynamic responses to brief periods of termed the hemodynamic response. Friston and col-
neural activity are delayed and dispersed in time. The leagues proposed in 1994 that an estimated hemo-
specific shape of these responses is of some importance dynamic response (treated as the impulse response
to the design and analysis of blood oxygenation level- function of a linear system) can be used to obtain a
dependent (BOLD), functional magnetic resonance im- predicted fMRI signal response for any arbitrary pattern
aging (fMRI) experiments. Using fMRI scanning, we of neural activity. These predicted signal responses can
examine here the characteristics and variability of then be used to test hypotheses regarding the effect of
hemodynamic responses from the central sulcus in
experimental treatments upon neural activity. Be-
human subjects during an event-related, simple reac-
cause this approach more accurately predicts the shape
tion time task. Specifically, we determine the contribu-
tion of subject, day, and scanning session (within a
of the fMRI signal (as compared to, for example, simply
day) to variability in the shape of evoked hemody- shifting a model of neural activity forward in time;
namic response. We find that while there is significant Bandettini et al., 1993), it affords greater statistical
and substantial variability in the shape of responses sensitivity and validity.
collected across subjects, responses collected during Based upon theoretical considerations, Friston and
multiple scans within a single subject are less variable. colleagues (1994) suggested that the shape of the
The results are discussed in terms of the impact of hemodynamic response can be modeled with a Poisson
response variability upon sensitivity and specificity of function. Later, Boynton and colleagues (1996) found
analyses of event-related fMRI designs. r 1998 Academic that a gamma function with two free parameters (plus
Press a pure phase delay) well modeled hemodynamic re-
sponses empirically derived from the primary visual
cortex of two subjects. Several groups (Courtney et al.,
1997; Dale and Buckner, 1997; Clark et al., 1998) now
An early and consistent observation regarding blood use the gamma model, and the parameter fits reported
oxygen level-dependent (BOLD), functional magnetic by Boynton and colleagues, to generate fMRI signal
resonance imaging (fMRI) data is that a sudden change predictions.
in neural activity produces a signal change that takes Because a single estimate of the hemodynamic re-
several seconds to develop and decay (Bandettini et al., sponse is used to analyze the data from different
1993). The sluggish nature of the BOLD fMRI signal is subjects, this approach assumes that any variability
a consequence of its hemodynamic origins: changes in that exists between subjects in hemodynamic response
neural activity engender changes in the local vascula- is minor. If this assumption is not true, then the general
ture and thus the local deoxyhemoglobin concentration, approach of using a ‘‘standard’’ hemodynamic response
to which BOLD is sensitive (Malonek and Grinvald, to analyze BOLD fMRI data from different subjects will
1996). Thus, BOLD fMRI provides a measure of the result in suboptimal power and perhaps invalid infer-
local, temporal pattern of neural activity, but only after ence. Several groups have now anecdotally noted that
that pattern has passed through a hemodynamic filter observed hemodynamic responses seem to vary from
that smoothes and delays the signal. Because the vast subject to subject (Boynton et al., 1996; Richter et al.,
majority of BOLD fMRI experiments test hypotheses 1996; Kim et al., 1997; Zarahn et al., 1997b). The
regarding neural activity, as opposed to vascular physi- purpose of the present report was to formally test the
ology, methods have been developed to account for the hypothesis that there exists variability in hemody-
temporal blurring imposed by the endogenous hemody- namic responses collected from different subjects. In
namic filter. addition, we evaluated other sources that might contrib-
The particular time-course of fMRI signal change ute to variability seen between subjects; namely, that

1053-8119/98 $25.00 360

Copyright r 1998 by Academic Press
All rights of reproduction in any form reserved.
 BOLD HEMODYNAMIC RESPONSES 361

within a subject hemodynamic responses might vary ment routine was used (part of SPM96b package;
from one scan to the next or from one day to the next. Friston et al., 1995b) without correction for ‘‘spin-
Measuring the relative contribution of these sources of history’’ (Friston et al., 1996). Next, a slice-wise motion
variability in hemodynamic response can inform as to compensation method was utilized that removed spa-
the physiological basis of response variability. tially coherent signal changes via the application of a
To test these hypotheses, a hemodynamic response partial correlation method to each slice in time (Zarahn
was obtained from each of several subjects. The sub- et al., 1997a).
jects performed a simple, event-related reaction time
task in which they made a bilateral button press every fMRI Datasets
16 s during fMRI scanning. This task was assumed to
produce a brief burst of neural activity within the Simple reaction time task datasets were obtained in
sensorimotor strip every 16 s, and the average fMRI a total of 41 (25 male) young [mean age (6SD) 5 23 6 3],
signal change that ensued after the button press events right-handed subjects. Subjects viewed a backlit projec-
was taken as an estimate of the hemodynamic response tion screen from within the magnet bore through a
for that subject. These responses were tested for the mirror mounted on the head coil. A white fixation cross
presence of significant variability. Additional subjects was constantly illuminated in the center of a black
participated in multiple scans, either on the same or background. Every 16 s the cross would briefly (500 ms)
different days. The responses obtained across scans for change to a white circle. The subject was instructed to
these subjects were also tested for variability to deter- monitor for this change, and to make a bilateral button
mine if the hemodynamic response is stable within a press with both thumbs on a fiberoptic game pad. A
subject across days or scans. total of 20 such trials were presented to each subject
during the scan.
A 16-s intertrial interval between button press events
METHODS
was selected as previous reports of BOLD fMRI re-
sponses indicate that the signal change has generally
MRI Technique and Initial Data Processing
run its course (i.e., has returned to baseline) after 16 s
Imaging was carried out on a 1.5T SIGNA scanner (e.g., Dale and Buckner, 1997). While residual signal
(GE Medical Systems) equipped with a fast gradient changes may linger after 16 s (e.g., Boynton et al.,
system for echoplanar imaging. A standard radiofre- 1996), these have been reported to be rather small in
quency (RF) head coil was used with foam padding to magnitude and would not be expected to greatly change
comfortably restrict head motion. High resolution sagit- the shape of the average response obtained here.
tal T1-weighted images were obtained in every subject. Nonetheless, the possibility that ‘‘overlap’’ of responses
A gradient-echo, echoplanar sequence was used to from one trial to the next might slightly alter the shape
acquire data sensitive to the BOLD signal at a TR 5 of the measured average response renders the hemody-
2000 ms, TE 5 50 ms. Resolution was 3.75 3 3.75 mm namic responses reported here suspect as ideal esti-
in plane, and 5 mm through plane, with no skip in between mates of the impulse response function of the system. It
planes (16 or 18 axial slices acquired). A total of 160 is important to note, however, that the possibility of
gradient-echo echoplanar images in time were obtained response overlap will not weaken or bias in any way the
per slice in each 320-s run. Twenty seconds of gradient and tests of variability of responses that are the focus of this
RF pulses preceded the actual data acquisition to allow report.
tissue to reach steady-state magnetization. Because of the regular spacing of the trials, sub-
Off-line data processing was performed on SUN jects were able to anticipate the occurrence of sti-
Sparc workstations using programs written in Interac- muli. Such anticipation might produce increases in
tive Data Language (Research Systems, Boulder, CO). neural activity within the motor cortex prior to the
After image reconstruction and prior to motion correc- onset of the stimulus (Georgopoulos et al., 1989).
tion, the data were sinc interpolated (by shifting the This possibility also makes it difficult to treat the
phase of the Fourier components) in time to correct for hemodynamic responses recorded here as good esti-
the interleaved fMRI acquisition sequence. This latter mates of the impulse response function, as the pre-
step is of particular importance here as hemodynamic cise moment of onset of neural activity cannot be
responses were to be compared across slices that were specified.
obtained at different points in the acquisition sequence Thirty-two of the subjects studied participated in
(and therefore at different points in time). If left only a single scan. One estimate of the hemodynamic
uncorrected, this would have introduced considerable response was obtained from each of these subjects.
variability and bias (a phase advance) into the hemody- Four other subjects participated in a total of five scans;
namic responses. The data were then motion corrected. each collected on a different day, spread out over a
First, a six parameter, rigid-body, least squares realign- several (5 to 11)-month period. From each of these
362 AGUIRRE, ZARAHN, AND D’ESPOSITO

subjects was obtained five estimates of the hemody- et al., 1997b) demonstrated voxel-wise false-positive
namic response. Finally, five additional subjects partici- rates in agreement with tabular values (data not
pated in five scans all collected during a single scanning shown).
session. One of these subjects did not possess signifi-
cant activity within the central sulcus during any of Derivation of Hemodynamic Responses
these scans. As a result, data from this subject could not
be used in the subsequent analyses and this subject is The central sulcus was defined upon each subject’s
not further discussed. From each of the remaining four T1 images by one of the authors (GA). The central
subjects, five estimates of the hemodynamic response sulcus was identified as the first medial-lateral sulcus
were obtained. posterior to, and not in contact with, the posterior
extent of the superior frontal sulcus on the superior
Creation of Statistical Maps most slices. The search volume included both the sulcus
and the surrounding gray matter, yielding a total (left
Voxel-wise analysis of the functional imaging data and right combined) search volume of ,200 voxels. The
was conducted to identify voxels with a significant statistical maps were thresholded at an F value corre-
response to the button-press events. Statistical maps sponding to a Bonferroni corrected, region-wise a 5
were created within the modified general linear model 0.05, and any voxels within the search region that
of Worsley and Friston (1995) using a Fourier basis set surpassed this threshold were identified. The average
of three sines and three cosines [frequencies (Hz) 5 time-series from this collection of suprathreshold vox-
0.0625, 0.125, 0.1875] (Josephs et al., 1997). (These six els was obtained and (i) filtered to remove frequencies
covariates provided a complete basis set for the eight below that of the paradigm and those around the
time points that they modeled as the data had been Nyquist (.0.244 Hz) and (ii) adjusted to remove the
filtered to remove the Nyquist frequency, and nuisance effects of nuisance covariates (Friston et al., 1995a).
covariates modeled the trial mean; see below.) Partial F The time series was then trial averaged and adjusted to
tests were used to evaluate the significance of the set the value of the first point to zero for display
variance in the data explained by these six covariates purposes. The resulting response was taken to be an
together. A specific advantage of this analysis approach estimate of the hemodynamic response of the system
is that sensitivity is not dependent on the shape of the for that scan.
response. Because all voxels with a significant response within
Previous work has shown that fMRI data are tempo- the central sulcus were averaged together, the studies
rally autocorrelated under the null-hypothesis (Aguirre conducted here are insensitive to variability in hemody-
et al., 1997a; Zarahn et al., 1997a). To account for this, a namic response from voxel to voxel within a region.
mean power spectrum was obtained from each dataset Such variability might be considered to arise from
by averaging the power spectra from each voxel. A differences in the diameter of the local vasculature
1/frequency (1/f ) function (Zarahn et al., 1997a) was within a voxel or the proximity of the voxel to neurally
then fit to this curve, ignoring those frequencies at active tissue (Menon et al., 1995). Future studies might
which power attributable to task might be expected. seek to address whether some subset of activated
The resulting curve was taken as an estimate of the voxels display more stable responses across scans
power spectrum of the data under the null-hypothesis. within a subject than other voxels. Nonetheless, the
The assumption underlying this analysis approach, possibility that voxel-by-voxel variability exists in the
that the noise is independent of stimulus temporal shape of the hemodynamic response within a subject
period, has been validated for fMRI (Boynton et al., will not impact the current studies that examine the
1996). variability of responses from the central sulcus region
The time-domain representation of the 1/f curve was as a whole.
placed within the K matrix (Worsley and Friston, 1995)
along with a filter designed to remove low frequency Statistical Analysis of IRF Variability
confounds (below 0.025 Hz) and high frequency noise at
and around the Nyquist frequency (above 0.244 Hz). It We considered three possible types of variability in
should be noted that these filtering components have evoked hemodynamic responses: variability (i) across
no effect upon the shape of the responses obtained, as subjects, (ii) within subject across days, and (iii) within
they impact frequencies that are either below that of subject and day across scans. To measure each type of
the task or above that passed by the hemodynamic variability, we obtained four sets of five hemodynamic
response of the system given a train of impulses as responses from each category:
input. Responses across subjects. Of the 32 hemodynamic
Application of this analysis approach to human BOLD responses obtained from different subjects, 20 were
fMRI data collected under the null-hypothesis (similar selected at random and divided into four groups of five.
to tests conducted in Aguirre et al., 1997a, 1998; Zarahn This provided for four tests of the variability of this
 BOLD HEMODYNAMIC RESPONSES 363

population of responses. Because these responses were A graphical example of the design matrices correspond-
collected from different subjects, one might attribute ing to a reduced and a full model is provided in Fig. 1.
variability in the shape of the response to the effect of The results of each hypothesis test were assessed by
subject. However, this dataset would also reflect any performing an appropriate partial F test comparing the
variability that exists between responses within the full and reduced models (Kirk, 1982, p. 179). Validation
same subject across days and any variability within of the model was performed using null-hypothesis
subject during a single scanning session. simulations in which five identical hemodynamic re-
Responses across days within subject. Four subjects sponses (plus normally distributed, computer gener-
were each scanned five times, each scan taking place on ated noise) served as the dependent data. These tests
a different day. Each of these sets of five responses confirmed that the test yields tabular false-positive
could be used to test for variability in responses ob- rates under the null-hypothesis.
tained from the same subject on different days. While The results of these tests of variability were 12 F
variability in response from one day to the next might values, 4 from each category (i.e., across subjects,
explain any variability observed in the responses, within subject across days, and within subject and
variability of responses within a scanning session and session across scans). An additional analysis deter-
subject would also be expected to contribute. mined if significantly greater variability (i.e., greater F
Responses across scans within subject and day. He- values) was observed in one group versus another.
modynamic responses were obtained from four subjects Because of the nonnormal distribution of the depen-
who each participated in five scanning runs during a dent data (i.e., F values) nonparametric tests were
single session on a single day. The source of variability employed. All 12 values were entered into a Kruskal–
in these data is scan-to-scan variability during a single Wallis analysis to determine if a significant effect of
category was present within the data. Then, pair-wise
session in a single subject.
Mann–Whitney tests were used to determine if greater
Thus, four sets of five hemodynamic responses were
F values were present within one given category as
obtained for each of three categories and tested for the
compared to another.
presence of significant variability of response. The
general approach adopted to test for variability was to
use partial F tests to compare regression models: a full
model that explicitly represented variability across
hemodynamic responses with interaction terms and a
reduced model that only modeled main effects of hemo-
dynamic response. Note that in all these analyses, the
hemodynamic response data (i.e., the dependent vari-
ables) were amplitude normalized to the sinc-interpo-
lated maximum positive excursion. This was done as
the hypothesized variability of interest was in the
shape of the hemodynamic response (relevant for its
use as a covariate) and not in its amplitude.
We wished to create a parsimonious set of indepen-
dent variables for use in the aforementioned tests of
variability. The remaining 12 of the 32 hemodynamic
responses obtained from different subjects were exam-
ined with a principal components analysis. The first
three eigenvectors of the principal components analysis
were then used as covariates within regression models
to test the hypotheses regarding sources of hemody-
namic response variability. The hemodynamic re-
sponses from each set were concatenated and served as
the dependent data. The reduced model contained FIG. 1. Example design matrices used to test for hemodynamic
response variability. (a) The reduced model. Arranged in columns
three main effect covariates, one for each eigenvector, from left to right, the model contains an overall intercept, three main
and trial effect covariates to mean center the residuals effect covariates, and four trial-effect covariates. Each main effect
of each response. The full model was identical to the covariate is composed of multiple, shifted versions of one of the
reduced model except for the addition of covariates eigenvectors derived from the central sulcus hemodynamic responses
which modeled interactions between the hypothesized from 12 subjects (see Fig. 3 and Table 1). (b) The full model. This
model is identical to that shown in a, except for the addition of
variability source and the behavior of the hemody- interaction terms designed to model the variability across hemody-
namic response. In these interaction covariates, each of namic responses. These two models were compared with a partial F
the three chosen eigenvectors was modeled separately. test (see Methods).
364 AGUIRRE, ZARAHN, AND D’ESPOSITO

RESULTS

IRF Descriptive Properties

Voxels that evidenced significant signal changes in
response to the button-press task were found within
the central sulci of all 32 of the subjects scanned once
[mean number of voxels (6SD) within central sulcus
region across subjects: 31 6 20, range: 6–91]. Across
subjects, the mean (6SD) peak signal change from
baseline was 2.05% 6 0.53 (n 5 32, range: 1.0 to 3.1%).
It is not possible to attribute a particular source to this
variability in response amplitude as it could be due to
differences in the amplitude of underlying neural activ-
ity (perhaps related to differences in amplitude of
motor response), differences in the transformation prop-
erties of the system, variations in nonphysiological
properties of the MRI scanner, or any combination of
these. As we are interested here in variability in the
filtering properties of the hemodynamic responses (i.e.,
the shape of the response), absent differences in raw
amplitude, all hemodynamic responses analyzed and FIG. 3. First three eigenvectors from a principle components
presented here have been normalized to the maximum analysis of central sulcus hemodynamic responses from 12 different
positive excursion. subjects. All vectors are orthonormal. The values of the eigenvectors
and their corresponding eigenvalues are provided in Table 1.
Twenty of the 32 hemodynamic responses obtained
from different subjects displayed some degree of de-
layed undershoot (i.e., signal values fell below baseline
postpeak), although there was some variability as to its shoot (Friston et al., 1994; Boynton et al., 1996) will be
degree [mean undershoot as a proportion of maximum unable to completely represent these empirically de-
positive excursion (6SD): 20.29 6 0.22 (n 5 20; range: rived responses.
20.02 to 20.81)]. It thus seems that models of the The sinc-interpolated time-to-peak was obtained for
hemodynamic response that do not include an under- all 32 hemodynamic responses. The mean (6SD) time-
to-peak of 4.7 6 1.1 s (n 5 32, range: 2.7 to 6.2) is in
rough agreement with that observed in previous stud-
ies of the BOLD hemodynamic system (Boynton et al.,
1996; Richter et al., 1996; Kim et al., 1997). Behavioral
(reaction time) data were obtained for 30 of the 32
subjects [mean RT (6SD) 5 368 6 68 ms]. There was
no discernible relationship between a subject’s time-to-
peak and reaction time [r (28 df ) 5 0.055, NS] (see Fig.
2) in agreement with previous studies (Richter et al.,
1996; Kim et al., 1997). While we might expect a linear,
unit slope relationship between these measures, our
failure to observe a significant correlation is unsurpris-
ing given the much larger variance in hemodynamic
time-to-peak relative to reaction time variance.

Generation of Eigenvectors of Evoked Hemodynamic

Responses
Twelve of the hemodynamic responses obtained from
subjects scanned a single time were selected at random.
These responses were subjected to a principle compo-
nents analysis. Figure 3 presents the first three eigen-
vectors from that analysis, which together explained
FIG. 2. Mean reaction time for each of 30 subjects versus the time
to peak of the hemodynamic response obtained from that subject’s
.98% of the variance within the set. The specific values
central sulcus. No significant relationship between the measures is of the three vectors and their eigenvalues are provided
present. in Table 1. As might be expected, the first principle
 BOLD HEMODYNAMIC RESPONSES 365

component strongly resembles the shape of hemody- TABLE 2

namic responses previously reported by several groups Results of Partial F Tests of Variability amongst
(Boynton et al., 1996; Richter et al., 1996). These Hemodynamic Responses
eigenvectors are used below in the analysis of variabil-
ity of hemodynamic responses. We also note that these Partial F tests (12, 20 df )
three vectors may serve as a parsimonious set of
Across Across Across
covariates in the analysis of other event-related fMRI
subject day scans
experiments.
Set 1 17.6 2.3 1.7
Tests of IRF Variability Set 2 20.2 9.6 0.6
Set 3 7.8 5.7 2.5
Several groups have anecdotally noted that hemody- Set 4 26.9 1.9 0.4
namic responses seem to differ from subject to subject Median 18.9 4 1.15
(e.g., Boynton et al., 1996; Kim et al., 1997) and that
responses appear to be more stable during a single Note. Four tests were conducted for each category of variability.
scanning session with a single subject (Kim et al., Significant (P , 0.05) responses in bold.
1997). We wished to test these observations and to
measure the variability observed in responses collected
within subjects on a single day, within subjects across group. The full model, but not the reduced model, also
days, and across subjects. contained interaction terms. These covariates modeled,
We asked, first, whether significant variability exists for each response, aspects of the response that were not
among hemodynamic responses collected from different shared across the group of five. The statistical test
subjects. Twenty hemodynamic responses obtained from measured if a significantly greater proportion of vari-
20 different subjects were randomly divided into four ance was explained by these interaction covariates
groups of five. The purpose of this division was to than would be expected by chance alone.
permit multiple tests of variability upon independent A partial F test revealed that the subject-by-response
datasets, thus allowing for replication. Each group of interaction covariates explained a significant amount
responses was tested for the presence of significant of variance for all four sets of responses (see Table 2).
variability. Variability was detected using partial F Thus, there are significant differences between hemody-
tests to compare statistical models that either explicitly namic responses collected across subjects from the
modeled the presence of variability among the re- same anatomical region. This test does not reveal,
sponses or did not. Both models contained three covari- however, whether or not this variability can be attrib-
ates (generated using the three eigenvectors) that uted to the effect of varying subject, day, and/or scan
modeled the ‘‘main effect’’ of responses across subjects session (as these are all confounded). Figure 4a presents
(see Fig. 1a). In effect, these covariates modeled the the set of across-subjects responses with the least
mean response, or the components of the hemodynamic variability while Fig. 4b presents the set with the
responses that were shared across subjects within the greatest variability.
Responses were also collected across multiple days
TABLE 1 from single subjects. In total, four subjects were each
scanned five times, with each scan occurring on a
Values of the Eigenvectors Shown in Fig. 2 different day over the course of several months. If
significant differences exist between these responses in
Principal components
a single subject, then scanning on different days can
Time (s) First Second Third contribute to response variability. Three of the four
subjects were found to have responses that varied
Eigenvectors significantly from one day to another. However, and
0 20.316 20.217 0.017 again, day was confounded here with ‘‘scan’’ as these
2 0.025 20.641 0.498 responses were obtained during different scanning
4 0.547 20.331 20.316 sessions. Figures 4c and 4d present two of these sets of
6 0.551 0.317 20.211
8 0.104 0.480 0.629
responses.
10 20.222 0.311 0.032 Finally, five hemodynamic responses were collected
12 20.335 0.067 20.331 from each of four subjects during a single scanning
14 20.353 0.013 20.319 session on a single day. The existence of significant
Eigenvalues differences between the responses obtained from a
0.743 0.222 0.018
single subject in this setting would suggest that hemo-
dynamic response variability exists from scan to scan.
Note. The eigenvalues have been normalized to the total variance. Partial F tests found significant variability in only one
366 AGUIRRE, ZARAHN, AND D’ESPOSITO

FIG. 4. Hemodynamic responses from the central sulcus, each sinc-interpolated and normalized to the maximum, interpolated positive
excursion. (a, b) Responses obtained from two different groups of five subjects each, representing the groups with (respectively) the least and
greatest variability in response shape. (c, d) Responses (least and greatest variability) collected from each of two subjects across different days.
(e, f) Responses (least and greatest variability) collected from each of two subjects across multiple scans in a single day. The F values reported
for each set of responses are the results of tests of variability (see Methods and Results). The key indicates the order of acquisition of the
responses for those collected within subjects (c–f ).
 BOLD HEMODYNAMIC RESPONSES 367

of the four subjects studied. Figures 4e and 4f present

two of these sets of responses from within subject and
day.
As can be seen from the median F scores (Table 2) for
each category of responses studied, the variability in
responses collected on different days is roughly 3.5
times greater than the variability present in responses
collected during a single scanning session. Moreover,
the median variability observed across subjects is over
16 times greater than that observed within a scanning
session. A Kruskal–Wallis test performed upon the F
values collected confirmed that there was a significant
effect of category (i.e., across subjects, across days,
across scans) upon the variability of responses observed
(P 5 0.02). Subsequent pair-wise tests revealed that
responses across subjects were significantly more vari-
able than either responses within a scanning session
(Mann–Whitney rank test, P 5 0.02) or responses across
days (Mann–Whitney rank test, P 5 0.04). Responses
across days, however, were not found to be significantly
more variable than responses within a scanning ses- FIG. 5. Three different summary models of the hemodynamic
response function. The thick solid line (circles) is the vector described
sion (Mann–Whitney rank test, P 5 0.08). Finally, it
by the gamma function of Boynton and colleagues (1996) [param-
should be noted that differences in variability of mean eters: t 5 1.25, n 5 3, d 5 2.5]. The thin solid line (diamonds) is a
RT across subjects compared to within subject cannot Poisson function [parameter 5 8 s] deconvolved with a Gaussian, s 5
account for these results, given the absence of any 1.9 s, to account for processing performed by Friston and colleagues
correlation between RT and time-to-peak observed (1994). The dotted line (squares) is the first eigenvector from Fig. 3.
above.
values, indicating a better fit between the model and
Response Variability and Standard Models of the
the population of responses. T tests performed upon the
Hemodynamic Response
t-transformed correlation values were used to assess
The significant variability in the shape of the hemody- the significance of these differences.
namic response function across subjects calls into ques- The Poisson model had a mean (6SD) correlation
tion the common use of a single, representative re- with the empirically obtained responses of 0.490 6
sponse function to model fMRI data from all subjects. 0.390 (n 5 20; range: 20.166 to 0.917). This indicates
Here we examined how well standard response func- that the Poisson function, on average, explained only
tions fit a population of responses from different sub- 25% of the variance present in the evoked responses. In
jects. We then tested if an alternative approach might contrast, the gamma model, and the first eigenvector
give rise to greater sensitivity: Given the stability of the reported here, both explained nearly 70% of the vari-
hemodynamic response observed within a single sub- ance in the evoked response on average. The average
ject, a promising analysis approach is to empirically (6SD) correlation of the gamma model with the popula-
derive a hemodynamic response for every subject, and tion of responses was 0.826 6 0.147 (n 5 20; range 5
then use that subject-specific response to analyze fur- 0.509 to 0.976) and the average for the first eigenvector
ther BOLD fMRI data from that subject. was 0.833 6 0.152 (n 5 20; range 5 0.501 to 0.991).
We examined three models of the hemodynamic The correlation values produced by both of these mod-
response (see Fig. 5). The first two of these are routinely els were significantly higher than those of the Poisson
used for the generation of covariates for the analysis of function [t(19 df ) 5 4.54, 3.87, P , 0.001].
BOLD fMRI data: the Poisson function described by Included in our datasets are four subjects who each
Friston and colleagues (1994) and the gamma function were scanned five times during a single scanning
described by Boynton and colleagues (1996). The third session, producing five hemodynamic responses per
model examined was the first eigenvector presented subject. We correlated the first hemodynamic response
here (Fig. 3). We measured the correlation between from each subject with the subsequent responses from
each of these models and each of the 20 hemodynamic that subject, resulting in four correlation values per
responses obtained from different subjects (recall that subject. The mean correlation value from each subject
the first eigenvector was generated from an indepen- was then taken to provide a set of four representative
dent data set). Of interest was the degree to which one correlation values, one from each subject, that describe
model or another tended to produce higher correlation the correspondence between a subject-specific hemody-
368 AGUIRRE, ZARAHN, AND D’ESPOSITO

namic response and subsequent responses from that perfectly valid, in that they will be unable to completely
subject during a scanning session. The mean (6SD) of model experimentally introduced variance. The particu-
these correlation values was 0.964 6 0.12 (n 5 4; lar impact of this error, as well as its degree, can be
range 5 0.947 to 0.976), indicating that using a subject- expected to vary greatly from one experimental design
specific response function results in a model that can to another (Aguirre and D’Esposito, 1998).
explain, on average, 92% of the variance in subsequent For traditional ‘‘blocked’’ fMRI designs, in which
evoked responses. The correlation values obtained from relatively long duration periods of neural activity are
this approach were significantly higher than those evoked, failure to perfectly model the shape of the
obtained from even the best of the standard models hemodynamic response will likely result in only small
tested just above [t(22 df ) 5 2.27, P 5 0.03]. The addi- decrements in sensitivity. Nonetheless, the use of a
tional explanatory power of this approach is not only more accurate response function can result in measur-
significant, but substantial as well—an additional 22% able increases in sensitivity, as has been demonstrated
of the variance of the hemodynamic response, on aver- for analyses of blocked fMRI experiments that used
age, was explained by using a subject-specific hemody- either a Poisson function or an empirically derived
namic response. estimate of the hemodynamic response (Aguirre et al.,
1997a).
For ‘‘event-related’’ fMRI designs, accurate estima-
DISCUSSION
tion of the hemodynamic response becomes more impor-
Given the stability of the hemodynamic response tant. For designs in which different trial types are
within subject relative to that across subjects, it seems randomly ordered (or perfectly counterbalanced) (e.g.,
that the particular shape of the hemodynamic response Dale and Buckner, 1997; Clark et al., 1998), failure to
within the central sulcal region is largely a conse- accurately specify the shape of the hemodynamic re-
quence of intersubject variability in physiology. The sponse can be expected to result in decrements in
smaller variability in responses within subject across sensitivity. This loss of power will likely increase as the
days suggests that this response may change over time spacing between the trials decreases (as the task
for some subjects. However, because subjects were not variance moves to ever higher frequencies) and might
studied during different scanning sessions on the same be substantial.
day (i.e., removed from the scanner and then reposi- It should be noted, however, that inaccurate charac-
tioned), it is not possible to unambiguously attribute terization of the hemodynamic response will not lead to
the variability in response across days to the effect of invalid inference under these circumstances. This is
day per se, as it is conceivable that different scanning because the random ordering of the trials makes any
sessions alone might introduce variability. difference between the assumed and actual hemody-
As has been noted, several groups make use of a namic response unbiased. This is to be contrasted with
priori estimates of the hemodynamic response of the event-related designs in which the different events of
BOLD fMRI system to test functional neuroimaging interest cannot be randomly ordered—classically, work-
hypotheses. This type of analysis approach can be ing-memory paradigms in which a delay period always
contrasted with alternative designs that use a flexible follows the presentation of a stimulus to be remem-
basis set instead. These approaches, like the Fourier bered (Courtney et al., 1997; Zarahn et al., 1997b). For
basis set design used here (Josephs et al., 1997), are designs of this kind, it is essential that the signal
notable for their equivalent sensitivity to any consis- change produced by one event type (e.g., the stimulus
tent pattern of signal change associated with neural presentation), be well modeled and not erroneously
events of interest. While the reduced assumptions attributed to an effect of another event type (e.g., the
inherent in such an analysis approach are sometimes delay period). Inferential false-positive results would
desirable, the lack of a model for the hemodynamic be expected if, for example, the actual hemodynamic
transform also limits the nature of inferences that can response of the subject was broader in time than
be drawn. This is because ‘‘significant’’ signal changes predicted by the assumed hemodynamic response. The
might be due to any one of numerous differences distribution of time-to-peak values in responses ob-
between two experimental conditions. Thus, a priori tained across subjects here (range: 2.7 to 6.2 s) suggests
estimates of the hemodynamic response will continue that this type of error of inference is a real possibility
to be necessary when focused hypotheses regarding the and should be guarded against by additional control
timing or intensity of neural activity are to be tested. conditions and tests in event-related designs of this
The presence of response variability across subjects kind (Zarahn et al., 1997b).
calls into question the use of a single representation of These considerations lead us to advocate the use of
a hemodynamic response to model evoked BOLD fMRI subject-specific hemodynamic responses for the analy-
signal changes in different subjects. Specifically, covari- sis of BOLD fMRI data. There is, however, an impor-
ates generated using a standard model will not be tant caveat to this suggestion. All of the responses that
 BOLD HEMODYNAMIC RESPONSES 369

were studied here were obtained from within the Boynton, G., Engel, S., Glover, G., Heeger, D. 1996. Linear systems
central sulcus. Anecdotal evidence suggests that the analysis of functional magnetic resonance imaging in human Vl. J.
Neurosci. 16:4207–4221.
shape of the hemodynamic response varies from one
Clark, V. P., Maisog, J. M., Haxby, J. V. 1998. An fMRI study of face
cortical region to another, although proof of this asser-
perception and memory using random stimulus sequences. J.
tion remains elusive, as it is unknown whether changes Neurophysiol., 79:3257–65.
in neural activity or vascular responsiveness underlie Courtney, S. M., Ungerleider, L. G., Keil, K., Haxby, J. V. 1997.
the regional response differences. Thus, the application Transient and sustained activity in a distributed neural system for
of a hemodynamic response derived from the primary human working memory. Nature 386:608–611.
motor cortex to test hypotheses in other areas of the Dale, A. M., Buckner, R. L. 1997. Selective averaging of rapidly
brain may also be strictly invalid. Future studies will presented individual trials using fMRI. Hum. Brain Map. 5:329–
be necessary to determine if (i) the hemodynamic 340.
response actually does vary from region to region and Friston, K., Ashburner, J., Frith, C., Poline, J.-B., Heather, J.,
(ii) if so, if this variability from region to region is Frackowiak, R. 1995b. Spatial registration and normalization of
images. Hum. Brain Map. 2:165–189.
greater than the variability within a region across
Friston, K., Williams, S., Howard, R., Frackowiak, R., Turner, R.
subjects. Alternatively, one might attempt to obtain an
1996. Movement-related effects in fMRI time-series. Magn. Reson.
estimate of the hemodynamic response from each corti- Med. 35:346–355.
cal area under study during preliminary studies. Friston, K. J., Holmes, A. P., Poline, J.-B., Grasby, P. J., Williams, S.
Finally, the stability of the response observed here C. R., Frackowiak, R. S. J., Turner, R. 1995a. Analysis of fMRI
within subject suggests the feasibility of fMRI designs time-series revisited. NeuroImage 2:45–53.
that attempt to detect small neural onset asynchronies. Friston, K. J., Jezzard, P., Turner, R. 1994. Analysis of functional MRI
Because the variability in the time-to-peak of hemody- time-series. Hum. Brain Map. 1:153–171.
namic responses collected from a single subject on a Georgopoulos, A. P., Crutcher, M. D., Schwartz, A. B. 1989. Cognitive
single day is rather minor, small (e.g., 50–100 ms) spatial motor processes: 3. Motor cortical prediction of movement
differences between two event types in the onset of direction during an instructed delay period. Exp. Brain Res.
75:183–194.
neural activity in a single region may be detected over
the course of repeated trials. Josephs, O., Rees, G., Turner, R., Friston, K. J. 1997. Event-related
fMRI. Hum. Brain Map. 5:243–248.
Kim, S.-G., Richter, W., Ugurbil, K. 1997. Limitations of temporal
ACKNOWLEDGMENTS resolution in functional MRI. Magn. Reson. Med. 37:631–636.
Kirk, R. E. 1982. Experimental Design: Procedures for the Behavioral
This work was supported by grants from the National Institutes of
Sciences. Brooks/Cole, Monterey, CA.
Health (NS01762 and AG13483) and the Charles A. Dana Foundation
and the American Federation for Aging Research. We also thank Malonek, D., and Grinvald, A. 1996. Interactions between electrical
Jessica Lease and Rajiv Singh for their help in data collection. activity and cortical microcirculation revealed by imaging spectros-
copy: Implications for functional brain mapping. Science 272:551–
554.
REFERENCES Menon, R. S., Ogawa, S., Hu, X., Strupp, J. P., Anderson, P., and
Ugurbil, K. 1995. BOLD based functional MRI at 4 tesla includes a
Aguirre, G. K., and D’Esposito, M. 1998. Experimental Design for
capillary bed contribution: Echo-planar imaging correlates with
Brain fMRI. In Functional MRI (P. A. Bandettini and C. Moonen,
previous optical imaging using intrinsic signals. Magn. Reson.
Eds.), in press. Springer Verlag, Berlin.
Med. 33:453–459.
Aguirre, G. K., Zarahn, E., and D’Esposito, M. 1998. A critique of the
Richter, W., Ugurbil, K., and Kim, S.-G. 1996. Limitations of temporal
use of the Kolmogorov-Smirnov (KS) statistic for the analysis of
resolution in fMRI. NeuroImage 3:S38.
BOLD fMRI data. Magn. Reson. Med. 39:500–505.
Aguirre, G. K., Zarahn, E., and D’Esposito, M. 1997a. Empirical Worsley, K. J., and Friston, K. J. 1995. Analysis of fMRI time-series
analyses of BOLD fMRI statistics. II. Spatially smoothed data revisited—Again. Neuroimage 2:173–182.
collected under null-hypothesis and experimental conditions. Neu- Zarahn, E., Aguirre, G. K., and D’Esposito, M. 1997b. A trial-based
roimage 5:199–212. experimental design for fMRI. NeuroImage 6:122–138.
Bandettini, P. A., Jesmanowitz, A., Wong, E. C., Hyde, J. S. 1993. Zarahn, E., Aguirre, G. K., and D’Esposito, M. 1997a. Empirical
Processing strategies for time-course data sets in functional MRI of analyses of BOLD fMRI statistics. I. Spatially unsmoothed data
the human brain. Magn. Reson. Med. 30:161–173. collected under null-hypothesis conditions. Neuroimage 5:179–197.