BOLD responses to stimuli: Dependence on frequency, stimulus form, amplitude, and repetition rate

www.elsevier.com/locate/ynimg NeuroImage 31 (2006) 585 – 599

BOLD responses to stimuli: Dependence on frequency, stimulus form, amplitude, and repetition rate P.A. Robinson, a,b P.M. Drysdale,a,b,* H. Van der Merwe, a E. Kyriakou, a M.K. Rigozzi, a,b B. Germanoska, a and C.J. Rennie a,b a

School of Physics, University of Sydney, NSW, 2006, Australia Brain Dynamics Center, Westmead Millenium Institute, Westmead Hospital and Western Clinical School of the University of Sydney, Westmead, NSW, 2145, Australia

b

Received 10 March 2005; revised 1 December 2005; accepted 20 December 2005 Available online 8 February 2006

A quantitative theory is developed for the relationship between stimulus and the resulting blood oxygen level-dependent (BOLD) functional MRI signal. The relationship of stimuli to neuronal activity during evoked responses is inferred from recent physiology-based quantitative modeling of evoked response potentials (ERPs). A hemodynamic model is then used to calculate the BOLD response to neuronal activity having the form of an impulse, a sinusoid, or an ERPlike damped sinusoid. Using the resulting equations, the BOLD response is analyzed for different forms, frequencies, and amplitudes of stimuli, in contrast with previous research, which has mostly concentrated on sustained stimuli. The BOLD frequency response is found to be closely linear in the parameter ranges of interest, with the form of a low-pass filter with a weak resonance at approximately 0.07 Hz. An improved BOLD impulse response is systematically obtained which includes initial dip and post-stimulus undershoot for some parameter ranges. It is found that the BOLD response depends strongly on the precise temporal course of the evoked neuronal activity, not just its peak value or typical amplitude. Indeed, for short stimuli, the linear BOLD response is closely proportional to the time-integrated activity change evoked by the stimulus, regardless of amplitude. It is concluded that there can be widely differing proportionalities between BOLD and peak activity, that this is the likely reason for the low level of correspondence seen experimentally between ERP sources and BOLD measurements and that non-BOLD measurements, such as ERPs, can be used to correct for this effect to obtain improved activity estimates. Finally, stimulus sequences that optimize the signal-to-noise ratio in event-related BOLD fMRI (efMRI) experiments are derived using the hemodynamic transfer function. D 2006 Elsevier Inc. All rights reserved.

* Corresponding author. E-mail address: [email protected] (P.M. Drysdale). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter D 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2005.12.026

Introduction The BOLD response is a measure of brain metabolic rate and, hence, is an indirect probe of neuronal activity. Current BOLD studies usually assume, either implicitly or explicitly, that (1) the relationship between neuronal spiking activity and BOLD is linear, and (2) that the coefficient of proportionality is always the same, at least for the same brain region. This underlies the common supposition that an increase in BOLD corresponds to an increase in stimulus. The first of these assumptions has only been partially explored to date (Friston et al., 2003; Janz et al., 2001), and the second does not appear to have been addressed in the literature. The BOLD response is usually measured in experimental protocols that employ a sustained stimulus that spans many seconds. However, recent experiments have begun investigating the event-related BOLD response (efMRI) in response to short stimuli in order to probe cognitive processing [see for example Miezin et al., 2000]. Also, it has been noted that the level of correspondence between ERP sources and BOLD or positron emission tomography (PET) measurements is remarkably low (Vitacco et al., 2002; Gamma et al., 2004). Hence, better understanding of the relationship between evoked neuronal responses and BOLD is necessary to extract accurate information on neuronal activity from the BOLD signal. Evoked response potentials (ERPs) are the brain’s transient electrical responses to discrete stimuli and are recorded from the scalp or from within the brain itself. These potentials are the product of extracellular current induced by neuronal activity (Niedermeyer and Lopes da Silva, 2002). Measurements of ERPs are commonly used in cognitive studies, providing high temporal resolution measurements of neuronal activity in the neocortex (Niedermeyer and Lopes da Silva, 2002). The higher time resolution of ERP experiments suggests they could be complementary to efMRI experiments which have good spatial resolution. Nonetheless, the connection between ERPs and BOLD remains ill defined. They measure different physiological

586

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

processes and thus could have differing utility for the purpose of inferring neuronal activity. This study explores these connections by combining an improved study of the activity – BOLD connection with recent insights into stimulus-activity-ERP links (Rennie et al., 2002). Recent physiologically based studies have modeled the BOLD as a response to neuronal activity. The hemodynamic study by Friston et al. (2003) provided a realistic approximation of the relationship between neuronal activity, vasodilatory response, blood inflow, and deoxyhemoglobin content, derived from the neuronal activity. While improving and correcting the Buxton et al. (1998) model, Obata et al. (2004) provided a refinement to the BOLD response of the Friston et al. (2003) model. The resulting Friston – Obata model extended the hemodynamic modeling of the so called Balloon Buxton et al. (1998) and Windkessel (Frank, 1899) models. Nunez (1981, 1995) and many others have studied the relationship between brain and scalp potentials in detail, taking into account geometry and volume conduction effects in the intervening tissues. More recently, physiologically based quantitative modeling has successfully reproduced ERPs, starting from basic physiological and anatomical data (Rennie et al., 2002; Robinson et al., 2003). The physiological and anatomical data for these quantitative studies are realistic and in agreement with existing literature for these parameters (Rowe et al., 2004; Robinson et al., 2004). Together, the results have established, in particular, a linear relationship between ERPs and neuronal activity. Hence, the ERP can be used to estimate and infer the form of the input to the Friston et al. (2003) equations assuming activity is proportional to metabolism, as evidenced by the work of Logothetis (2002). It should be stressed that this result differs strongly from some early suggestions that ERP peaks and troughs result from high excitatory and inhibitory activity, respectively—the two activities vary very nearly in phase (Rennie et al., 2002). The present study investigates the relationship between ERPs and BOLD response using the ERPs to infer the event-related neuronal input to the hemodynamic state equations of Friston et al. (2003) and Obata et al. (2004), which are then used to calculate the BOLD response. This paper thus quantitatively interrelates the various aspects of neuronal activity, including the time course of activity, and the resulting ERP and BOLD response. The Friston – Obata hemodynamic model used here employs physiologically based hemodynamics to link stimulus and neuronal activity, to BOLD response. This paper analyzes a fundamental nonlinear physics model in contrast to Dynamic Causal Modelling (Friston et al., 2003) a bilinear statistical model of hemodynamics. The model presented here complements Statistical Parametric Mapping (SPM) analyses, which are based on empirical descriptions of the BOLD response [e.g., Friston, 2004]. Our approach also enables a more fundamental and unified understanding of a variety of phenomena and of why the hemodynamic response function has the empirically observed form. In Section 2, the hemodynamic and ERP aspects of the model are introduced. Section 3 studies the response of the BOLD signal to variation in the frequency of sinusoidally varying neuronal activity. Section 4 explores the dependence of the BOLD response on variations in the form and amplitude of a stimulus that results in an ERP. Section 5 analyzes the response to short pulses and derives the most appropriate choice of stimulus patterns to optimize efMRI response.

Theory The relationship between ERPs and BOLD response can be modeled theoretically using a compact set of differential equations. This study combines the (Friston et al., 2003) model that links the BOLD response to neuronal activity with key insights from a mathematical model of that activity and makes use of the resulting equations to describe the BOLD response and how it relates to ERPs. Hemodynamic state equations Friston et al. (2000) described the hemodynamic responses, as derived from neuronal activity using the Balloon/Windkessel model (Buxton et al., 1998; Frank, 1899). Their equations connect the neuronal activity z, to changes in the vasodilatory response s, blood inflow F, change in blood volume m, and deoxyhemoglobin content q. Each location (e.g., a voxel) is considered independent of all others, which is a reasonable approximation so long as they are not too small. The equations for a given voxel are ds ¼ z  js  cð F  1Þ; dt

ð1Þ

dF ¼ s; dt

ð2Þ

s

dr ¼ F  m1=a ; dt

ð3Þ

s

dq FEð F;qÞ m1=a q ¼  : dt q m

ð4Þ

where a is Grubb’s exponent (Grubb et al., 1974), s is the hemodynamic transit time, c is the flow-dependent elimination constant, j is the rate of signal decay, and E is expressed as E ð F;qÞ ¼ 1  ð1  qÞ1=F ;

ð5Þ

where q is the resting oxygen extraction fraction (Buxton et al., 1998). Eq. (1) has been corrected as it had inconsistent dimensionality in Friston et al. (2003), whereas in Friston et al. (2000), s had dimensionality of inverse time (s1), and F was nondimensional. To preserve consistency with these definitions, the flow-dependent elimination constant c is defined here to have the dimensionality of inverse squared time (s2), rather than being a rate with units of (s1) as stated in Friston et al. (2003). The BOLD signal that results from the above dynamics is denoted by y and is a function of the blood volume and deoxyhemoglobin content. Obata et al. (2004) corrected an error and slightly redefined parameters in the BOLD response equation originally appearing in Buxton et al. (1998). The Friston et al. (2003) paper appeared before the correction was published, so it used the older form of the response function given by y ¼ V0 ½k1F ð1  qÞ þ k2F ð1  q=mÞ þ k3F ð1  mÞ;

ð6Þ

where k 1F, k 2F, and k 3F are earlier BOLD response parameters defined in Buxton et al. (1998) and depend on the magnetic field

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

and the resting oxygen extraction fraction q. Table 1 shows estimates of these parameters for a 1.5 T field. In this paper, we use the corrected form of the BOLD response given by Obata et al. (2004) y ¼ V0 ½k1 ð1  qÞ þ k2 ðm  qÞ  k3 ð1  mÞ;

ð7Þ

where k 1, k 2, and k 3 are dimensionless parameters depending on the magnetic field and the resting oxygen extraction fraction q. For a 1.5 T field, these parameters have been estimated by Obata et al. (2004) and are given in Table 1. Table 1 also shows the other parameter values Friston et al. (2003) used in their modeling. Basic ERP and physiology The previous section does not take into account the biological underpinning of the neuronal activity. However, recent physiologically based models have successfully reproduced ERPs (Rennie et al., 2002). One key outcome is that a linear relationship between the ERP and neuronal activity exists over a wide range. Hence, the ERP can be used to estimate and infer the form and time course of neuronal activity. This estimated neuronal input can then be substituted into the Friston et al. (2003) equations and used to predict corresponding efMRI responses. The theoretical and practical issues of defining what constitutes neuronal activity in BOLD models were discussed by Glaser et al. (2004). The model presented here relies on two inferences. First, that there is a linear relationship between neuronal activity expressed via the local field potential and metabolism and thus blood flow. This inference is supported by the experimental work of Logothetis (2002) and is also used in the BOLD model of Friston et al. (2000). Second, that there is a linear relationship between ERPs and underlying neuronal activity, based on the work of Rennie et al. (2002), for example. In contrast with previous analyses where the neuronal input has usually been modeled by a ‘‘boxcar’’ function that represents a sustained stimulus, in this paper, the neuronal input is inferred from typical forms of the ERPs observed for discrete stimuli. Measured ERPs exhibit a wide range of signal forms but a number of general properties can be stated. For short duration discrete stimuli, ERPs consist of Fearly_ components (20 – 60 ms after the stimulus) which are the product of initial processing of the stimulus and the Flate_ components which begin about 100 ms after the onset of the stimulus, resemble EEGs in amplitude (10 AV), have typical frequency content Table 1 Parameter values used by Friston et al. (2003), plus the quantities k 1 – k 3 introduced by Obata et al. (2004) Parameter

Description

Value

j c s a q V0 k1 k2 k3 k 1F k 2F k 3F

Rate of signal decay Flow-dependent elimination constant Hemodynamic transit time Grubb’s exponent Resting oxygen extraction fraction Resting blood volume fraction Obata et al.’s BOLD response coefficients

0.65 s1 0.41 s2 0.98 s 0.32 0.34 0.02 7q 1.43q 0.43 7q 2 2q – 0.2

Friston et al.’s BOLD response coefficients

587

2 – 40 Hz and consist of oscillations which die away within 500 ms. ERPs can thus be approximated by damped sinusoids. Indeed, the use of damped sinusoids in ERP analysis has a long history in decomposition of ERPs into components and phenomenological modeling (Franaszczuk and Blinowska, 1985; Freeman, 1964). In Section 4, a damped sinusoid ERP approximation is used for simplicity to investigate the hemodynamic response for a variety of input characteristics. The mathematical form used is zðt Þ ¼ aebt sinðbdtÞ;

ð8Þ

for t > 0, and z(t) = 0 for t < 0, where a is the amplitude, the positive constant b is the rate of decay (the inverse of the damping time), and the positive quantity D is proportional to the number of oscillations per damping time. If bY0, but d increases so that bd is a finite, nonzero quantity, a semi-infinite sine wave is produced, suitable for studying the relationship of steady-state evoked potentials (SSEP) to the corresponding BOLD response. We stress again the critical result that excitatory and inhibitory activities vary almost in phase during ERPs, so the time course of the ERP closely reflects that of the total neuronal activity (Rennie et al., 2002).

Frequency dependence of the BOLD response The first aspect of BOLD response to be considered is its dependence on the frequency of the neuronal input z(t). This can be obtained by evaluating the transfer function of the linearized hemodynamic response equations. The limitations of the linear analysis are considered in Section 4. Theory Eqs. (1) – (4) are linearized by first noting that for steady-state conditions, the equations must satisfy ds/dt = 0, dF/dt = 0, dm/dt = 0, and dq/dt = 0. These yield steady-state values for z (0) = 0, of s (0) = 0, F (0) = 1, m (0) = 1, q (0) = 1, and y (0) = 0. Linearizing about these steady-state values yields dsð1Þ ¼ zð1Þ  jsð1Þ  cF ð1Þ ; dt

ð9Þ

dFð1Þ ¼ sð1Þ ; dt

ð10Þ

drð1Þ 1 1 ¼ F ð1Þ  mð1Þ ; s as dt

ð11Þ

dqð1Þ q þ ð1  qÞlnð1  qÞ ð1Þ 1  a ð1Þ 1 ð1Þ F  m  q ; ¼ qs as s dt

ð12Þ

h   i yð1Þ ¼ V0  k1 qð1Þ þ k2 mð1Þ  qð1Þ þ k3 mð1Þ ;

ð13Þ

where all quantities with superscripts are linear perturbations relative to the steady-state values. We now determine the transfer function y (1)(x)/z (1)x) that relates the BOLD response to neuronal activity at a given angular frequency x, with x = 2pf, where f is the usual frequency in Hertz.

588

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

This is evaluated by obtaining the transfer functions relating each of the state variables s (1), F (1), m (1), and q (1) to the neuronal input z (1). We first note that Eqs. (9) and (10) only involve s (1) and F (1), while (11) and (12) are both driven by F (1). Eqs. (9) and (10) correspond to a damped harmonic oscillator in F (1) (Friston et al., 2000) and yield transfer functions Tsz ðxÞ ¼

TFz ðxÞ ¼

sð1Þ ðxÞ ix ; ¼ 2 zðxÞ x þ ijx  c F

shown the optimal 14-s spacing between stimulations observed in efMRI experiments (Bandettini and Cox, 2000) is simply the reciprocal of this resonant frequency. At higher frequencies, the transfer functions given by Eqs. (14) and (15) fall off as f 1 and f 2. The transfer functions relating m (1) and q (1) to blood inflow F (1) can be obtained from Eqs. (11) and (12), giving

ð14Þ TmF ðxÞ ¼

mð1Þ ðxÞ a ; ¼ 1  isax F ð1Þ ðxÞ

TqF ðxÞ ¼

qð1Þ ðxÞ q þ ð1  qÞlnð1  qÞ ¼ qð1  isxÞ F ð1Þ ðxÞ

ð1Þ

ðx Þ 1 : ¼ 2 zðxÞ x þ ijx  c

ð17Þ

ð15Þ

For simplicity, transfer functions are denoted as T xy = x (1) ( f)/y (1) ( f) in this paper. Figs. 1a and b show )Tsz ( f))2 and )T Fz ( f))2, respectively. Each transfer function shows a resonant peak, with the peak in F (1) corresponding to the vasodilatory system’s resonant frequency  1=2 fp ¼ c  j2 =2 =2p: ð16Þ For the parameters in Table 1, this corresponds to f p = 0.07 Hz, giving a period of approximately 14 s. In Section 5.3, it is



1a : ð1  isxÞð1  isaxÞ

ð18Þ

These transfer functions are seen in Figs. 1c and d, showing peak responses at f = 0. The transfer function T rF is a relatively wideband low pass filter of the blood inflow response, whereas T qF is a narrower bandwidth low pass filter.

Fig. 1. Frequency dependence of squared transfer functions relating state variables and neuronal inputs for parameters in Table 1. (a) )Tsz )2, (b) )T Fz )2, (c) )T rF )2, (d) )T qF )2, (e) )T rz )2, (f) )T qz )2.

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

589

The transfer functions for m (1) and q (1) can now be expressed in terms of the neuronal activity: Tmz ðxÞ ¼

mð1Þ ðxÞ a ; ¼ 2 ðx þ ijx  cÞð1  isaxÞ zðx Þ

ð19Þ

Tqz ðxÞ ¼

qð1Þ ðxÞ ða þ b  1  iabsxÞ ; ¼ zðx Þ ð1  isxÞð1  isaxÞðx2 þ ijx  cÞ

ð20Þ

where b¼

q þ ð1  qÞlnð1  qÞ : q

ð21Þ

Fig. 1e shows the transfer function relating m (1) to the neuronal activity. As might be expected from the weakness of the filtering shown in Fig. 1c, the response of m (1) to neuronal activity in Fig. 1e is very similar to the blood inflow response in Fig. 1b. The response approximately peaks at the natural frequency of the vasodilatory response f p . However, at higher frequencies, the response falls off as f -3, rather than the slower f -2 fall off of the blood inflow response. Fig. 1f shows the response of deoxyhemoglobin content q (1) to neuronal activity. The low pass filtering effect shown in Fig. 1d reduces the strength of the resonance in the transfer function of q (1) when compared to the blood inflow response, and the peak response frequency is shifted slightly lower. From Eq. (7), the BOLD response is seen to depend on changes in blood volume and deoxyhemoglobin content. Using Eqs. (7), (19), and (20), we thus find the overall transfer functions that relate the BOLD signal to the flow rate and neuronal activity: TyF ðxÞ ¼ V0

T yz ðxÞ¼ V0

aðk2 þ k3 Þð1  isxÞ  ðk1 þ k2 Þða þ b  1  isabxÞ ; ð1  isxÞð1  isaxÞ ð22Þ aðk2 þ k3 Þðisx  1Þ þ ðk1 þ k2 Þða þ b  1  iabsxÞ ; ð1  isxÞð1  isaxÞðx2 þ ijx  cÞ ð23Þ

where b is given by Eq. (21). The transfer function (23) relating BOLD response to neuronal activity is shown in Fig. 2. The properties of the BOLD response can be deduced from the transfer function (23) and from Fig. 2. The BOLD signal responds relatively effectively to low-frequency variations in the neuronal input z (1)(t) but not to high-frequency changes, leading to the f 6 decline in the squared transfer function )Tyz )2 at higher frequencies. The transfer function has a weak resonance at f p  0.07 Hz for the parameters in Table 1. The dependence of the BOLD transfer function on its physiological parameters is now considered. The position of the resonance, and the sharp decline of the transfer function depends strongly only on the rates of signal decay j and flow-dependent elimination c, while the response at very low frequencies f  0 is sensitive to the resting oxygen extraction fraction q. For the typical range 0.2 < q < 0.55 (Friston et al., 2000), the zero-frequency response remains fairly constant relative to the peak BOLD frequency response. Friston et al. (2000) suggested very high q > 0.7 may be present in the visual cortex. For q in this very high range, the response near f  0 falls, slightly sharpening the peak of the transfer function at 0.07 Hz, as seen in Fig. 2. It should be noted that the uncorrected form of the BOLD response, Eq. (6), has a

Fig. 2. Transfer function relating the BOLD response to neuronal activity )Tyz (x))2 as given by Eq. (23) for the parameters in Table 1 and two other values of q. Resting oxygen extraction fractions are q = 0.34 as per Table 1 (solid line), q = 0.7 (dashed line) and q = 0.22 (dotted line).

significantly different frequency response at high q from that described here, with greater sharpening of the resonance at high q. Experimentally, it has been observed that Grubb’s exponent varies between its steady-state value and its value under stimulation. With 6 s of stimulation, a falls from a  0.32 to a  0.18 (Mandeville et al., 1999). The shape of the transfer function is relatively insensitive to this variation but the magnitude of the transfer function increases by approximately 30% at all frequencies. Comparison with other studies Only a few experimental studies of the BOLD frequency response are available in the literature. For example, Nunez and Silberstein (2000) obtained an approximate frequency response by inverting an analytical empirical fit of the experimental BOLD impulse response. Their approximate response was investigated for frequencies 1 – 20 Hz, well into the tail of the BOLD frequency response, beyond the primary frequency range of interest in this paper (approximately 0 – 0.5 Hz). The tail of the frequency response in Nunez and Silberstein (2000) varied as x 4, in contrast to x 3 obtained via the more detailed analysis in this paper. Logothetis (2002) presented (in his Fig. 19) BOLD responses to square wave stimuli at four different frequencies for a monkey brain in connection with in situ field measurements of neuronal activity. The square wave stimulus is dominated by the principal sine wave component, and the response to higher harmonics is severely attenuated by the hemodynamics in the frequency range investigated by Logothetis (2002). Hence, the data can be used to infer four points on the frequency response curve at frequencies corresponding to the stimulus repetition rates. To process the data for comparison with (23), the published curves were scanned and digitized, linear trends were removed, and the resulting BOLD response between 20 and 50 s after stimulus onset was Fourier transformed in each case. (This temporal window was selected to exclude the bulk of the transient associated with stimulus onset.) The resulting estimates of the transfer function have large uncertainties due to the short time interval covered by the data and their intrinsic noisiness. These points are shown in Fig. 3. The data are at significantly higher frequencies than the shoulder of the transfer function and noise contributes strongly to the observed signal.

590

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

Fig. 3. Comparison of theoretical and experimental BOLD frequency response. The four data points are experimental BOLD frequency response derived from Logothetis (2002) data. The solid line is the theoretical BOLD frequency response (23) with added noise derived from Logothetis (2002) data. The dotted line corresponds to theoretical BOLD frequency response without ambient noise.

Because noise contributes significantly to the data, the estimated ambient noise is added to the theoretical frequency response to allow more direct comparison with experiment. The ambient noise is estimated in the following manner: each of the four time series is Fourier transformed to yield a frequency spectrum, the noise components of the four resulting power spectra are averaged (i.e., neglecting the peak at the stimulation frequency in each case), and the resulting mean noise is fitted linearly vs. frequency to estimate the mean ambient noise spectrum. The theoretical transfer function with this additional estimated noise is also shown in Fig. 3. The final experimental and theoretical estimates in Fig. 3 are consistent, but the comparison is inconclusive because it is dominated by noise due to the frequencies available lying well above the pass band of the BOLD response (which has f | 0.15 Hz), which results in a very small signal. In future, data for frequencies f < 0.15 Hz will be required to determine the transfer function experimentally and make detailed comparisons with theory.

Response to model ERPs In this section, variation in the BOLD response due to changes in the form and amplitude of the neuronal signal z(t) are investigated for realistic, damped sinusoidal neuronal activities inferred from typical experimentally observed ERPs which are closely proportional to neuronal activities (Rennie et al., 2002). Unlike the previous frequency response analysis in Section 3, the analysis in this section considers the full nonlinear Friston et al. (2000) hemodynamic model with the Obata et al. (2004) correction, as given by Eqs. (1 – 4, rather than the linearized model embodied in Eqs. (9) – (12). This analysis also enables the nonlinear limits of the hemodynamic model to be probed. Before considering the nonlinear form of the equations, the linear system response obtained in the previous section can be used to make a number of predictions of the form and amplitude of the neuronal signal and the resulting BOLD response.

As discussed above, experimentally observed ERPs have major frequency components in the range 2 – 40 Hz. Since neuronal activity is closely proportional to the ERPs, it can be inferred that neuronal activities due to short duration stimuli will also have major frequency components in the range 2 – 40 Hz (Rennie et al., 2002). As was seen in Fig. 2, the linear BOLD response is a low pass filter with a shoulder around 0.1 Hz, effectively cutting out at least an order of magnitude below the dominant frequency content of neuronal input inferred from ERPs. We thus predict that the BOLD response depends only the (nearly) zero-frequency components of the neuronal input. In mathematical terms, the BOLD response in the linear model is thus approximately proportional to the signed area under the neuronal activity curve (i.e., the time integral of the change in neuronal activity due to the stimulus) for physiologically plausible short neuronal activity changes inferred from measured ERPs. The analysis of the BOLD response to a damped sinusoidal input is deferred to Section 5.2 where we confirm in detail the proportionality between BOLD responses and the signed area under the neuronal activity curve. The full nonlinear hemodynamic model embodied in Eqs. (1 – 4 is numerically integrated for cases where the driving input z(t) is given by a damped sinusoid (8). The relationships between BOLD response and the parameters a, b, and d of the damped sinusoid (8) are determined by varying each parameter independently. The BOLD response to very large stimuli is also explored to establish the limit beyond which linearity is not a good approximation. Fig. 4a shows the BOLD response as a function of amplitude a for fixed decay rate b and oscillation frequency bd. The amplitude of neuronal input is linearly related to the BOLD response, at least in the range 0  a  2 shown. This linear relationship is consistent with the nonlinear model being close to linear in the physiologically plausible regime and the experimentally observed BOLD response being approximately linear for short stimuli. In fact, as seen in Fig. 4b, the amplitude must increase to a  80 for the nonlinear model to show strong nonlinearity between increased amplitude and BOLD response for damped sinusoids. This is far beyond the physiologically reasonable parameter range of 0  a  1 discussed above. The situation changes substantially if the shape and form of the neuronal signal are changed by varying b and/or d. For example, Fig. 4c shows the BOLD response is proportional to 1/b. The signed area under the neuronal activity curve is also proportional to 1/b, confirming the BOLD response is proportional to the signed area under the neuronal activity curve regardless of changes in the decay constant b. Fig. 4d shows the BOLD response as a function of d. Fig. 4e shows the ratio of the BOLD response to the timeintegrated change in neuronal activity (the signed area under the activity curve) for varying d. As predicted, the BOLD response is again proportional to this signed area. Overall, in agreement with the linear analysis in the previous section, the BOLD response is found to be closely proportional to the signed area under the neuronal activity curve (i.e., proportional to the low-frequency components) for a in the physiologically realistic regime. This relationship applies to variations in the shape and form of the neuronal signal due to variation of b and d. The BOLD response is thus only proportional to the amplitude of the neuronal activity signal so long as shape and form of this signal do not vary. This qualification is of significant practical importance because it has usually been implicitly assumed that the relationship between neuronal spiking activity and BOLD is linear and the coefficient of proportionality is always the same, at least for the

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

591

Fig. 4. Maximum BOLD response of nonlinear model for damped sinusoid stimulation Eq. (8). Parameters are those given in Table 1 BOLD response is shown with variation in (a) physiological range of a, (b) extremely large a, (c) 1b, (d) d. (e) The Maximum BOLD response divided by the signed area of neuronal activity, showing that the ratio of these quantities is almost exactly constant.

same brain region. Our results show this is not true except under the highly restrictive condition that the neuronal activity only changes by a normalization factor, keeping a constant time profile in all other respects. (This is certainly not true in most cases where evoked responses to different stimuli are compared and can only be established by non-BOLD methods in any case.) The BOLD response for short stimuli is proportional to the signed area under the neuronal activity signal. Different fMRIs are not directly comparable with regard to strength of the response even at the same point in the brain except if strength of response is defined in the sense of two events having comparable signed area under neuronal activity curve. Fig. 5 illustrates the above key results. Two pairs of neuronal activity input and the resulting BOLD responses are shown, which were obtained by numerical integration of the nonlinear Friston et al. (2000) model with the Obata et al. (2004) correction. Although the peak amplitude of the neuronal activity is greater in Fig. 5c than a, the resulting BOLD response is half as large in panel (d) as in (b) reflecting the proportionality between BOLD response and

the time integrated neuronal activity. This result has not previously been suspected because previous analyses have considered only long stimuli and/or boxcar stimulus models which cannot reveal these features. By using short stimulus activities inferred from ERPs including ones with different properties and negative segments, this property is now evident. The practical consequences of the above result are as follows: ERPs are known to show a very wide range of shapes and forms and the forms change with different kinds of stimuli. Thus, the fact that two different stimuli yield very different, or even oppositely signed, BOLD responses in a particular brain region cannot be used to infer that the peak amplitude of the neuronal activity differs in proportion to the BOLD responses in the two cases. Indeed, it is not even possible to infer which stimulus has greater activity, except in a time integrated sense. This may be the underlying reason that correspondences between experimental ERP sources and BOLD or PET signals (the latter also having a long time constant) are remarkably weak (Gamma et al., 2004; Vitacco et al., 2002).

592

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

Fig. 5. Neuronal activities (a) and (c) and the corresponding BOLD response (b) and (d) respectively for the nonlinear model. Parameters are as given in Table 1.

The differences between the BOLD responses represent changes in the very low-frequency components (signed areas) between the different neuronal activity curves. In order to be able to determine directly whether a stronger BOLD response corresponds to greater neuronal activity in cognitive studies, some external non-BOLD measurements (e.g., ERPs) must be available to indicate that neuronal activities have not changed in shape and form. Alternatively, where the neuronal activities do differ in time profile, measurements of the ERP profiles can be used to determine the relative proportionalities between BOLD and activity in the various cases and thereby correct for this effect to extract improved neuronal activity information from the BOLD responses.

The integral (24) can be evaluated by contour integration and yields an impulse response of yð1Þ ðt Þ ¼ Aet=s  Bet=as þ ejt=2 ½Csinðx0 t Þ  Dcosðx0 tÞ; where x0 ¼ given by

ð1  bÞðk1 þ k2 Þs ; 1  js þ s2 c

ð26Þ

B ¼ V0

a2 ðk1  k3 Þs ; 1  ajs þ a2 s2 c

ð27Þ

C¼

Impulse response The general form of the impulse response y (1)(t) is obtained from the transfer function (23), via the inverse Fourier Transform yð1Þ ðtÞ ¼

1 2p

Z

V

V

yð1Þ ðxÞ ð1Þ z ðxÞeixt dx; zð1Þ ðxÞ

ð24Þ

where z (1)(x) = 1 for a delta function impulse input z (1)(t) = d(t).

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c  j2 =4 and the constants A, B, C, and D are

A ¼ V0

Impulse and ERP-related responses and optimized repetition patterns We next calculate the optimal temporal pattern with which to present short stimuli in order to increase the BOLD signal to noise ratio in cognitive studies, which is often low. In the previous section, the nonlinear hemodynamic model was numerically analyzed for neuronal activity input with physiologically plausible parameters inferred from ERPs. The results confirmed that the linear model is a good approximation of the nonlinear model over the physiologically realistic parameter range. Consequently, the linear model is used in this section.

ð25Þ

  1 ð2  ajsÞ ; Að2  jsÞ  B 2sx0 a

D ¼ A  B:

ð28Þ ð29Þ

For the parameters in Table 1, A = 0.060, B = 0.0047, C = 0.051, D = 0.055, x 0 = 0.55 s1, and the decay times are s = 0.98 s, as = 0.31 s, and 2/j = 3.1 s. The impulse response is dominated by the damped oscillatory (i.e., the third) term of (25) after ¨2 s since the first term has decayed by this time, and the second term of (25) is smaller to begin with and decays even faster. Fig. 6 shows the analytic impulse response curves for two different sets of parameters. The impulse response in Fig. 6 with the high resting oxygen extraction fraction q shows the initial dip (t < 1 s) observed in some studies (Jezzard et al., 1997; Logothetis, 2002). Friston et al. (2000) explored this dip by numerically integrating Eqs. (1) – (4) with the earlier form of the BOLD Eq. (6). With the updated form (7), the initial dip is still present, but its size is smaller for the same parameter values. The analytic form (25) confirms the initial dip is only strongly sensitive to q, k 1, k 2, and k 3. Since the k i are dependent on

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

Fig. 6. BOLD response y(t) under impulse neuronal stimulus as given by Eq. (25) for parameters in Table 1 except for resting oxygen extraction fraction q which is given by q = 0.34 (solid line) and q = 0.75 (dotted line).

magnetic field, the analytic impulse response shows increased initial dip with increasing field strength, as suggested by Obata et al. (2004) in their development of the updated BOLD Eq. (7). Thus, at higher magnetic fields, the BOLD response is more sensitive to deoxygenation prior to flow changes. The principal features of the analytic impulse response are also consistent with experimentally measured BOLD responses to short stimuli (Logothetis, 2002; Miezin et al., 2000). The short stimulus analog of the post-stimulus undershoot, known from the long stimulus literature (Jezzard et al., 2001; Kwong et al., 1992), is also seen after t = 7 s in Fig. 6; its duration and profile are consistent with long stimulus analogs. Many empirically fitted BOLD impulse responses appear in the literature. These range from individual gamma variate (Friston et al., 1994) and Poisson distributions (Bandettini and Cox, 2000) to fits based on sums of orthogonal basis sets including sets derived from Fourier sets (Josephs et al., 1997), ‘‘canonical’’ sets (Friston et al., 1998), and sets of gamma functions of different widths and lags (Boynton et al., 1996). The impulse response derived above differs from all these approaches by being parameterized in terms of underlying physiological parameters. To illustrate similarities and differences, we compare several widely used empirical equations with the impulse response derived above. For example, Friston et al. (1994) used a gamma-variate distribution. Although the form of the gamma-variate distribution stated in that paper was in a discrete-time form, it can be expressed in a continuous time form given by yð1Þ ðtÞ ¼

kt=t0 ek ; Cðt=t0 þ 1Þ

593

form (25) shows the initial dip and post-stimulus undershoot. The narrower pulse width of our form leads to the higher relative pulse amplitude of the analytic form compared to the empirical forms. The narrower pulse width is consistent with experimental data published by Miezin et al. (2000) in their Figs. 2 and 7). It should be noted both initial dip and post-stimulus undershoot can appear in impulse responses obtained using basis set techniques (Friston et al., 1998; Josephs et al., 1997). The impulse response (25) is of comparable mathematical simplicity to existing ‘‘single element’’ impulse responses (Bandettini and Cox, 2000; Friston et al., 1994) but has the advantage of being described in terms of the underlying physiological parameters, instead of being fitted empirically. Dynamics are thus related naturally to these physiological parameters, and we thus argue that this systematically derived form of the impulse response should be used in future modeling to achieve improved accuracy and insight. One area where the new impulse response will be valuable is to replace simple ‘‘single element’’ descriptions of the BOLD impulse response. It remains appropriate for SPM analyses to use basis set methods to closely fit empirical impulse response functions, but future work in this direction may benefit from the physiologically based insights obtained here into the appropriate form of the basis functions. Josephs et al. (1997) fitted the impulse response function using a basis set of Fourier terms using sixteen complex frequencies up to 0.5 Hz. Our hemodynamic transfer function provides why their approach was successful: frequency components up to 0.5 Hz are required to capture the high-frequency response of the impulse response function, while a frequency spacing of 0.03 Hz is required to resolve the dominant low-frequency part of the hemodynamic transfer function. Response to ERP-like damped sinusoidal input The analytic form of the BOLD response y (1)(t) to a damped sinusoidal ERP-like activity curve (8) can be derived analogously to the impulse response. In this case, z (1)(x) is given by zð1Þ ðxÞ ¼

abd ðb þ ixÞ2  b2 d2

;

ð32Þ

which is substituted into (24) to obtain y (1)(t). The result can be expressed as a sum of exponentially decaying sinusoids, but we do not show it explicitly here as the full result is complicated and

ð30Þ

where k = 7.69, t 0 = 1 s, and the Euler function C has been used to generalize the distribution to continuous time. The Euler function satisfies C(1 + n) = n!, for nonnegative integers n and generalizes this factorial function to real n (Kreyszig, 1993). Bandettini and Cox (2000) used a Poisson form given by yð1Þ ðtÞ ¼ ðt=t0 Þb et=c ;

ð31Þ

where b = 8.6, c = 0.55 s, and t 0 = 1 s. Fig. 7 compares the impulse response derived in Eq. (25) with the empirical approximations (30) and (31). For comparison purposes, the signed areas of the empirical approximations (30) and (31) are normalized to the signed area of the analytic response (25). All three profiles have a strong positive peak, but only our

Fig. 7. Impulse response derived from Friston – Obata hemodynamic model with parameters given in Table 1 (solid line) compared with empirically derived gamma-variate and Poisson process impulse responses (dashed and dotted lines respectively). All three models have been normalized to the same signed area for comparison purposes.

594

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

an excellent, much simpler, approximation can be easily obtained since the sinusoidal frequency bd/2p in (8) is much larger than 0.07 Hz in the cases of interest. The approximation can be systematically derived from (32) and (24) by noting that the integral that gives the analytic BOLD response can be written in terms of a contour integral. The value of a contour integral is, by Cauchy’s theorem, related to the sum of the residues at the poles of the integrand which is analytic (Kreyszig, 1993). In the present case, the integrand has two additional poles relative to the corresponding one for the impulse response. This results in two additional exponentially decaying sinusoids, but these can be neglected in the parameter range of interest because their residues turn out to be small relative to those of the other poles. The remaining residues are equal to those in the impulse-response case, but multiplied by a factor of approximately z (1)(x))x = 0 (the signed area under the neuronal activity curve). Thus, we conclude that the BOLD response due to the damped sinusoid neuronal activity is obtained simply by multiplying the (unit) impulse response (25) by the signed area under the neuronal activity curve. This result is not unexpected, given the low-pass characteristics of the BOLD transfer function, and is exemplified by the similarity between the BOLD responses in Fig. 5, which are in response to damped sinusoidal activity, and the BOLD impulse response in Fig. 6. Optimized repetition patterns Low signal to noise ratios (SNRs) are a significant problem in measuring BOLD responses to short stimuli. Consequently, use of an optimized sequence of stimuli is very important for obtaining an adequate BOLD response. The optimization of stimulus timing based on measured, empirically fitted, or idealized hemodynamic response functions has been extensively explored in the literature (Josephs and Henson, 1999; Friston et al., 1999; Liu et al., 2001; Buckner et al., 1998; Birn et al., 2002; Burrock et al., 1998). Broadly, such studies have concluded that blocked stimulus patterns are superior to other designs, especially from a detection power perspective. Further, these studies showed that different blocked designs had similarly strong BOLD responses (Liu et al., 2001). Given the importance of optimizing stimulus patterns to improve BOLD response, we now systematically derive an optimal temporal pattern of stimuli from the Friston – Obata hemodynamic model. By using a physiologically based model, we thus gain improved insight into optimization of stimulus patterns. It should be noted that an individual neuronal activity event inferred from an ERP is of short duration 0.5 s and thus has a frequency spectrum approximately an order of magnitude broader than the BOLD response transfer function, implying weak coupling to the BOLD response for a single short stimulus. This accounts for the significant signal to noise problems encountered for short stimuli. However, by using a sequence of such stimuli, the activity spectrum can be reshaped to be better concentrated in the optimal part of the transfer function. This enables the SNR to be improved, as we show next. Traditionally BOLD measurements are made with long stimulus durations (SDs) of order 14 – 20 s. In order to correct for experimental noise, these are interspersed with similarly long interstimulus intervals (ISIs) to ensure a 50% – on 50% – off stimulus pattern. By comparing with the frequency spectrum of such long pulses with the BOLD transfer function, it is found that

the zero-frequency and fundamental modes of the stimulus induced activity lie near the peak of the transfer function, thereby ensuring a strong BOLD response. The relatively uniform low-frequency response of the BOLD transfer function ensures that the BOLD response will be strong as long as the fundamental lies close to its peak (i.e., SD a 12 s). In addition, the BOLD response is somewhat stronger for SD a 36 s when the harmonic (at three times the fundamental frequency of the stimulus repetition rate) also contributes significantly to the BOLD response. Unfortunately, long stimulus durations are poorly suited to cognitive studies because cognitive processing times are short, so short stimuli are needed to probe them. Such stimuli are dominated by highfrequency content, which results in small BOLD response, and the SNR is then a problem. In order to improve the SNR, Bandettini and Cox (2000) experimentally and theoretically investigated optimal interstimulus intervals (ISI) for short stimulus durations (SDs). Using the empirical response (31), they found that BOLD response was maximal when SD + ISI 12 – 14 s, in agreement with their experimental results. The form of the BOLD transfer function obtained in Section 3 immediately demonstrates why this pattern is optimal. Our analysis shows for this ISI, one of the principal frequency components of the stimulation pattern lies at the weak resonance at 0.07 Hz, leading to the stronger BOLD response than for other ISIs. We now consider the more general problem of optimizing a general pattern of stimuli and ISIs to increase BOLD response. The purpose here is to exploit the hemodynamics to maximize the BOLD response, subject to the constraint that the individual stimulation events lead to independent neuronal responses. We increase the BOLD response by optimizing the temporal pattern of stimuli to concentrate neuronal activity at frequencies near the peak of the BOLD transfer function. Concentrating the neuronal activity near the peak of the BOLD transfer function is the optimal manner for reshaping the neuronal activity frequency spectrum to maximize BOLD response. The optimization problem will now be stated mathematically. The neuronal activity of a sequence of independent neuronal activity responses is given by X zð1Þ ðt Þ ¼ g ðt  tk Þ; ð33Þ k

where g(t) is the neuronal activity in response to a single stimulation event, k is the number of stimulation events and t k their onset times. Both k and the t k are as yet unknown. The choice of t k is constrained as each neuronal activity should remain independent. We assume that there exists a minimum time T required to ensure independence of neuronal responses, and thus, the minimum permissible time between stimuli; the constraint t k  t k-1  T is thus required. The frequency spectrum of the neuronal activity is given by the Fourier transform Z X F ðx Þ ¼ g ðt  tk Þeixt dt: ð34Þ k

Because of the existence of the BOLD resonance, the optimization problem reduces to finding the number and onset times t k of stimuli that maximize the magnitude )F(x m )) of the frequency spectrum at the resonant frequency x m . The frequency x m of the maximal BOLD transfer Tyz is slightly less than the frequency x p of the maximum of T Fz , obtainable from (16).

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

Eq. (34) can be simplified by first reversing the order of the sum and integral, which yields F ðxÞ ¼ ~ Xg ðt  tk Þeixt dt:

ð35Þ

k

Next, the integration variable in each of the integrals is changed to u k = t  t k , giving F ðxÞ ¼ ~ Xg ðuk Þeixðuk þtk Þ duk :

ð36Þ

k

The common integral term is then factored out of the sum and the subscript k dropped off the dummy variables u k without loss of generality. This yields F ðxÞ ¼ Xg ðuÞeixu du ~ eixtk ;

ð37Þ

k

¼ g ðxÞ ~ eixtk ;

ð38Þ

k

where g(x ) is the Fourier Transform of g(t). To maximize the BOLD response, the pattern of stimuli should be chosen to maximize the magnitude of the sum in (38) at x = x m . [The factor g(x m) is a fixed property of the individual neuronal response g(t) and thus cannot be varied in the optimization process.] We next deduce the optimal stimulus pattern in two steps: (i) the optimal pattern is shown to be periodic with period T m = 2p/x m (ii) the optimal pattern within each oscillation period T m is derived. Both results can be derived by considering the sum in Eq. (38) geometrically at x = x m . Each term corresponds to a unit vector in the complex plane and the onset times t k determine the phase terms of each of the unit vectors. The magnitude )F(x m )) of the sum is the magnitude of the sum of these vectors. Thus, the optimal stimulation pattern must have each stimulus in the same complex half plane to avoid the tendency for vectors to cancel each other. Since each unit vector has a phase term periodic with period T m = 2p/x m , this corresponds to optimizing the stimulus pattern within a sequence of half-periods of duration T m /2 separated by nonstimulation intervals of the same length. The periodicity of the optimal stimulus pattern can now be deduced. For short stimulus cognitive studies, ERPs have durations well under 1 s, so T  1 s is appropriate and thus T < T m /2. The assumption that T < T m /2 guarantees that the neuronal activity spectrum can be significantly reshaped. Noting that the only constraint on the times of stimuli is t k  t k1  T and since T < T m /2, the pattern of stimuli within a given half-period of stimulation is independent of the arrangement in the previous half-period containing stimuli. Since the stimulus pattern in each half-period can thus be optimized independently (and thus are all the same), the optimal pattern of stimuli is periodic with period T m . This is a stronger result than that described in the previous paragraph, which only showed that the optimal stimulus pattern occurs in T m /2 blocks. Here, each block of stimulus has been shown to be independent, and therefore, the optimal pattern is periodic. The optimal pattern of stimulation within one half-period is now deduced. It was noted above that to ensure independence of neuronal responses, the t k must be selected so that t k  t k1  T. This constraint is geometrically equivalent to imposing a minimum angle of x m T between successive vectors. Geometrically, )F(x m )) is maximized if the t k is chosen to satisfy the following two conditions: (i) all the vectors lie in the same half plane (i.e., within a half-period T m/2 of the optimal oscillation), as described above, and (ii) the largest possible number of vectors are placed within

595

this half plane. The latter condition is subject to the constraint imposing a minimal angle between the vectors which implies a maximum N stim = T m /2T number of vectors that can be placed within a half plane. Fig. 8 shows the optimal stimulation pattern schematically. Within the optimization literature, the experimental results of Bandettini and Cox (2000) are of interest here because we demonstrate that they show possible evidence of a measurable effect on the optimal stimulation patterns arising from the weak resonance in the hemodynamic transfer function. The Bandettini and Cox (2000) scheme stimulated approximately every 14 s. This corresponds to stimulating at the maximum of the BOLD transfer function. However, this scheme only stimulates at the phase where the system is maximally responsive; i.e., only at the peak of the sinusoid in Fig. 8. In contrast, the optimal scheme inserts as many additional stimuli as possible in the same half period subject to the constraint that they are at least a time T apart. The interval T is determined by the form of the particular stimulus and neuronal response and is typically ¨1 s in cognitive studies. The optimal scheme thus involves stimulating as often as possible for 7 s of every 14 s and repeating this pattern every 14 s. A simple analogy is that of pushing a swing: if one is restricted to delivering very short pushes in a single direction, maximum amplitude is achieved if these are delivered as often as possible during the half-cycles when the swing is moving away. There is nothing fundamental about any particular zero of time, but once swinging has commenced, the relevant half-cycles are unambiguously defined. This last point clarifies the nature of the schematic sinusoid in Fig. 8, whose phase is only defined once the hemodynamic response commences, but which then determines unambiguous half cycles in which stimuli should be placed to optimize SNR. Significantly, stimulation in both half cycles will actually reduce the response and worsen the SNR. The optimal nature of this stimulation pattern is directly demonstrated in Figs. 9 and 10. In both figures, the full nonlinear Friston – Obata model is integrated for the shown stimulation patterns, and the resulting BOLD signal is displayed. Fig. 9 shows that the optimal stimulation pattern period T m is 14 s. For T m < 14 s, the BOLD signal has reduced amplitude. For T m > 14 s, the BOLD signal is not optimal in the presence of noise, although the BOLD signal amplitude is similar. This is since when comparing signals with similar amplitude, a BOLD response which spends more equal times at high and low response is more resistant to the effects of noise and thus superior

Fig. 8. Schematic of optimal sequence of stimuli for maximizing BOLD response while keeping short duration neuronal activities independent. The sinusoidal curve schematically represents the hemodynamic system’s instantaneous responsivity to neuronal activity over one period T m . For the physiological parameters in Table 1, T m  14 s.

596

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

Fig. 9. Stimulation patterns (arrows), inferred neuronal activity (solid line) and resulting predicted BOLD response (dashed line). The predicted BOLD response is derived from integration of the full nonlinear Friston – Obata model. Neuronal activities due to individual stimulus events are modeled using a damped sinusoid Eq. (8) with a = 1, b = 10, and c = 2.5. The stimulation pattern period T m is (a) 7 s, (b) 14 s, and (c) 21 s.

to less balanced BOLD responses (see for example Liu et al. (2001), Birn et al. (2002)). For similar reasons, Fig. 10 shows that the optimal stimulation pattern fills the first half cycle of each period T m . The optimal stimulation pattern we have derived can now be considered within the context of existing optimization literature. The optimal stimulation pattern derivable from the physiologically based hemodynamic model is consistent with existing optimization studies which suggest blocked designs are superior. The weak hemodynamic resonance has not been previously noted in the literature. This resonance is predicted to select a particular blocked pattern as optimal (the 7 s stimulating 7 s off pattern described above). Bandettini and Cox’s experimental results relating to ISIs suggest the optimal pattern predicted above will show stronger BOLD response in experiments. Further experimental confirmation of the weak resonance is thus a priority.

Visualizing the BOLD response as a hemodynamic resonance provides a valuable method for understanding why the strength of BOLD response depends only really on moderate variations in the repeat frequency of the blocked design. For frequencies up to ¨0.08 Hz, the hemodynamic transfer function is relatively constant (see Fig. 2), so reshaping the neuronal activity spectrum to have dominant frequencies in this range will result in similar strength BOLD responses. Existing optimization studies often convolve an idealized impulse response function with stimulus patterns to determine optimal stimulus patterns. A variety of different functional forms are assumed for the response functions (e.g., gamma variate or Poisson distributions). Although in each case blocked designs are found to be superior, the choice of idealized response function can have consequences for which blocked patterns are found to be optimal. Further analysis of Eq. (24) shows that impulse response functions without post-stimulus

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

597

Fig. 10. Stimulation patterns (arrows), inferred neuronal activity (solid line) and resulting predicted BOLD response (dashed line). The predicted BOLD response is derived from integration of the full nonlinear Friston – Obata model. The number of stimulation events within a ¨14-s stimulation pattern is varied between frames. Neuronal activities due to individual stimulus events are modeled using a damped sinusoid Eq. (8) with a = 1, b = 10, and c = 2.5.

undershoot do not exhibit weak resonance, so inclusion of this undershoot will ensure that predicted optimizations are in closer agreement with experiment.

Summary and discussion We have combined insights from quantitative analysis of ERPs with hemodynamic modeling to infer the connection between stimuli and the BOLD response they evoke. The hemodynamic modeling is based on fundamental physiology, rather than empirical approaches commonly found in the literature. In particular, we have systematically explored the BOLD response evoked by impulsive, sinusoidal, and ERP-like damped-sinusoidal activities and have determined the pattern of short stimuli that maximizes the BOLD response. The results also

unify existing literature and provide insight into how the underlying physiological mechanisms result in the observed BOLD response and why blocked patterns are superior for optimizing BOLD response. For the first time, the BOLD transfer function is obtained, and the impulse response function is expressed in terms of the underlying physiological parameters of the hemodynamic response. The BOLD frequency response is found to act as a low-pass filter with a weak resonance at approximately 0.07 Hz for the relevant physiological parameters (see Table 1). The hemodynamics were found to be very nearly linear for physiologically plausible neuronal activities. The dependence of the linear BOLD response on the physiological parameters was then systematically derived from the model. The BOLD impulse response and BOLD response due to a damped sinusoidal stimulus were obtained analytically and found to be expressible as compact sums of exponentially decaying

598

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599

sinusoids. This analytic expression for the response exhibits the initial dip and post-stimulus undershoot, which are both absent in some commonly used empirical approximations to the impulse response; it can thus be expected to be more accurate in applications. Future work improving basis set functions for SPM style analyses may benefit from the physiologically based insights obtained here into the appropriate form of the basis set functions. A key result of this analysis is that the BOLD response to short stimuli depends crucially on the time course of the activity, not just on its peak value or typical amplitude. Indeed, for short stimuli, the BOLD response is proportional to the time-integrated evoked activity change relative to background, rather than any other features of the activity. This is likely to account for the low level of correspondence observed between ERP sources and BOLD or PET signals. The above dependence on the temporal course of the activity implies that, unless it is known that the profile of the neuronal activity is the same constant for different stimuli, the BOLD response alone cannot be used to infer their relative evoked neuronal activities, except in the above time-integrated sense. If peak neuronal activity is to be measured, additional non-BOLD information is required to determine the time profile of the activity and correct for this effect to obtain accurate activity estimates. ERP data, although obtained at larger spatial scales, are likely to be particularly useful in this context. The problem of BOLD response signal to noise ratios for short stimuli has been addressed here via our understanding of the BOLD transfer function as a weakly resonant low-pass filter with a resonant period of circa 14 s. The optimal stimulation pattern for short stimuli is to repeat them as many times as possible for the first half of each period, subject to requirement that neuronal activities evoked by successive stimuli be independent, with no stimuli in the second half of each period. This improves on the Bandettini and Cox (2000) scheme of stimulating every 14 s, which turns out also to exploit the hemodynamic resonance at ¨0.07 Hz, albeit to a lesser extent. Future empirical studies of optimization of BOLD response should appropriately model poststimulus undershoot as it has been found here the optimal pattern depends on presence or absence of this undershoot.

Acknowledgments We thank Mathew Barton and Leanne Williams for the stimulating discussions. The Australian Research Council supported this work.

References Bandettini, P.A., Cox, R.W., 2000. Event-related fMRI contrast when using constant interstimulus interval: theory and experiment. MRM 43, 540 – 548. Birn, R.M., Cox, R.W., Bandettini, P.A., 2002. Detection versus estimation in event-related fMRI: choosing the optimal stimulus timing. NeuroImage 15, 252 – 264. Boynton, G.M., Engel, S.A., Glover, G.H., Heeger, D.J., 1996. Linear systems analysis of functional magnetic resonance imaging in human V1. J. Neurosci. 16, 4207 – 4221. Buckner, R.L., Goodman, J., Burrock, M., Rotte, M., Koustaal, W., Schacter, D., Rosen, B., Dale, A.M., 1998. Functional – anatomic

correlates of object priming in humans revealed by rapid presentation event-related fMRI. Neuron 20, 285 – 296. Burrock, M.A., Buckner, R.L., Woldorff, M.G., Rosen, B.R., Dale, A.M., 1998. Randomized event-related experimental designs allow for extremely rapid presentation rates using functional MRI. NeuroReport 9, 3735 – 3739. Buxton, R.B., Wong, E.C., Frank, L.R., 1998. Dynamics of blood flow and oxygenation changes during brain activation: The Balloon model. MRM 39, 855 – 864. Franaszczuk, P.J., Blinowska, K.J., 1985. Linear model of brain electrical activity—EEG as a superposition of damped oscillatory modes. Biol. Cybern. 53, 19 – 25. Frank, O., 1899. Die Grundform des arteriellen Pulses. Z. Biol. 85, 91 – 130. Freeman, W.J., 1964. Use of digital adaptive filters for measuring pyriform evoked potentials from cats. Exp. Neurol. 10, 475 – 492. Friston, K.J., 2004. Experimental design and statistical parametric mapping. In: Frackowiak, R.S.J., et al., (Eds.), Human Brain Function, (2nd edR) Academic, Amsterdam, pp. 599 – 632. Friston, K.J., Jezzard, P., Turner, R., 1994. Analysis of functional MRI time-series. Hum. Brain Mapp. 1, 153 – 171. Friston, K.J., Fletcher, P., Josephs, O., Holmes, A., Rugg, M.D., Turner, R., 1998. Event-related fMRI: characterizing differential responses. NeuroImage 7, 30 – 40. Friston, K.J., Zarahn, E., Josephs, O., Henson, R.N.A., Dale, A.M., 1999. Stochastic designs in event-related fMRI. NeuroImage 10, 607 – 619. Friston, K.J., Mechelli, A., Turner, R., Price, C.J., 2000. Nonlinear responses in fMRI: The Balloon model. Voltera kernels and other hemodynamics. NeuroImage 12, 466 – 477. Friston, K.J., Harrison, L., Penny, W., 2003. Dynamic causal modelling. NeuroImage 19, 1273 – 1302. Gamma, A., Lehmann, D., Frei, E., Iwata, K., Pascual-Marqui, R., Vollenweider, D., 2004. Comparison of simultaneously recorded [H15 2 O]-PET and LORETA during cognitive and pharmacological activation 2004. Hum. Brain Mapp. 22, 83 – 96. Glaser, D., Friston, K.J., Mechelli, A., Turner, R., Price, C., 2004. Haemodynamic modelling. In: Frackowiak, R.S.J., et al., (Eds.), Human Brain Function, (2nd edR) Academic, Amsterdam, pp. 823 – 841. Grubb, R.L., Raichle, M.E., Eichling, J.O., Ter-Pogossian, M.M., 1974. The effects of changes in Pa CO2 on cerebral blood volume, blood flow, and vascular mean transit time. Stroke 5, 630 – 639. Janz, C., Heinrich, S.P., Bach, M., Hennig, J., 2001. Coupling of neural activity and BOLD fMRI response: new insights by combination of fMRI and VEP experiments in transition from single events to continuous stimulation. MRM 46, 482 – 486. Jezzard, P., Rauschecker, J.P., Malonek, D., 1997. An in vivo model for functional MRI in cat visual cortex. MRM 38, 699 – 705. Jezzard, P., Matthews, P.M., Smith, S.M. (Eds.), Functional MRI An Introduction to Methods. Oxford, Oxford, pp. 161 – 162. Josephs, O., Henson, R.N.A., 1999. Event-related functional magnetic resonance imaging: modelling, inference and optimization. Philos. Trans. R. Soc. London 354 (1215), 1228. Josephs, O., Turner, R., Friston, K.J., 1997. Event-related fMRI. Hum. Brain Mapp. 5, 243 – 248. Kreyszig, E., 1993. Advanced Engineering Mathematics, 7th edR Wiley, New York. Kwong, K.K., Belliveau, J.W., Chesler, D.A., Goldberg, I.E., Weisskoff, R.M., Poncelet, B.P., Kennedy, D.N., Hoppel, B.E., Cohen, M.S., Turner, R., Cheng, H.-.M., Brady, T.J., Rosen, B.R., 1992. Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc. Natl. Acad. Sci. U. S. A. 89, 5675 – 5679. Logothetis, N.K., 2002. The neural basis of the blood-oxygen-leveldependent function magnetic resonance imaging signal. Philos. Trans. R. Soc. Lond., B 357, 1003 – 1037. Liu, T.T., Frank, L.R., Wong, E.C., Buxton, R.B., 2001. Detection power, estimation efficiency, and predictability in event-related fMRI. NeuroImage 13, 759 – 773.

P.A. Robinson et al. / NeuroImage 31 (2006) 585 – 599 Mandeville, J.B., Marota, J.J.A., Ayata, C., Zaharchuk, G., Moskowitz, M.A., Rosen, B.R., Wiesskoff, R.M., 1999. Evidence of cerebrovascular postarteriole Windkessel with delayed compliance. J. Cereb. Blood Flow Metab. 19, 679 – 689. Miezin, F.M., Maccotta, J.M., Ollinger, J.M., Petersen, S.E., Buckner, R.L., 2000. Characterizing the hemodynamic response: effects of presentation rate, sampling procedure, and the possibility of ordering brain activity based on relative timing. NeuroImage 11, 735 – 759. Niedermeyer, E., Lopes da Silva, F., 2002. Electroencephalography: Basic principles, Clinical Applications and Related Fields, 4th edR Williams and Wilkins, Baltimore, USA. Nunez, P.L., 1981. Neocortical Dynamics and Human EEG Rhythms. Oxford, Oxford. Nunez, P.L., 1995. Neocortical Dynamics and Human EEG Rhythms. Oxford, Oxford. Nunez, P.L., Silberstein, R.B., 2000. On the relationship of synaptic activity to macroscopic measurements: does co-registration of EEG with fMRI make sense? Brain Topogr. 13, 79 – 96.

599

Obata, T., Liu, T., Miller, K.L., Luh, W., Wong, E.C., Frank, L.R., Buxton, R.B., 2004. Discrepancies between BOLD and flow dynamics in primary and supplementary motor areas: application of the balloon model to the interpretation of BOLD transients. NeuroImage 21, 144 – 153. Rennie, C.J., Robinson, P.A., Wright, J.J., 2002. Unified neurophysical model of EEG spectra and evoked potentials. Biol. Cybern. 86, 457 – 471. Robinson, P.A., Rennie, C.J., Rowe, D.L., O_Connor, S.C., Wright, J.J., Gordon, E., Whitehouse, R.W., 2003. Neurophysical modeling of brain dynamics. Neuropsychopharmacology 28, S74 – S79. Robinson, P.A., Rennie, C.J., Rowe, D.L., O_Connor, S.C., 2004. Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Hum. Brain Mapp. 23, 53 – 72. Rowe, D.L., Robinson, P.A., Rennie, C.J., 2004. Estimation of neurophysiological parameters from the waking EEG using a biophysical model of brain dynamics. J. Theor. Biol. 231, 413 – 433. Vitacco, D., Brandeis, D., Pascual-Marqui, R., Martin, E., 2002. Correspondence of event-related potential tomography and functional magnetic resonance imaging during language processing. Hum. Brain Mapp. 17, 4 – 12.