Comparison of hemodynamic response nonlinearity across primary cortical areas

www.elsevier.com/locate/ynimg NeuroImage 22 (2004) 1117 – 1127

Comparison of hemodynamic response nonlinearity across primary cortical areas David A. Soltysik, a,* Kyung K. Peck, b,1 Keith D. White, c Bruce Crosson, d and Richard W. Briggs b,2 a

Department of Nuclear and Radiological Engineering, University of Florida, Gainesville, FL 32610, USA Department of Radiology, University of Florida, Gainesville, FL 32610, USA c Department of Psychology, University of Florida, Gainesville, FL 32611, USA d Department of Clinical and Health Psychology, University of Florida, Gainesville, FL 32610, USA b

Received 21 October 2003; revised 3 March 2004; accepted 8 March 2004

Hemodynamic responses to auditory and visual stimuli and motor tasks were assessed for the nonlinearity of response in each of the respective primary cortices. Five stimulus or task durations were used (1, 2, 4, 8, and 16 s), and five male subjects (aged 19 F 1.9 years) were imaged. Two tests of linearity were conducted. The first test consisted of using BOLD responses to short stimuli to predict responses to longer stimuli. The second test consisted of fitting ideal impulse response functions to the observed responses for each event duration. Both methods show that the extent of the nonlinearity varies across cortices. Results for the second method indicate that the hemodynamic response is nonlinear for stimuli less than 10 s in the primary auditory cortex, nonlinear for tasks less than 7 s in the primary motor cortex, and nonlinear for stimuli less than 3 s in the primary visual cortex. In addition, neural adaptation functions were characterized that could model the observed nonlinearities. D 2004 Elsevier Inc. All rights reserved. Keywords: BOLD fMRI; Linearity; Stimulus duration; Auditory; Motor; Visual; Neural adaptation

Introduction Functional magnetic resonance imaging (fMRI) using blood oxygenation level-dependent (BOLD) contrast has emerged as a noninvasive method to measure vascular oxygenation changes due to brain activity (Kwong et al., 1992; Ogawa et al., 1992). (For a review, see Hennig et al., 2003). The shape of the evolving hemodynamic response is a function of stimulus intensity, stimulus * Corresponding author. Department of Biophysics, Medical College of Wisconsin, 8701 Watertown Plank Road, Milwaukee, WI 53226. E-mail address: [email protected] (D.A. Soltysik). 1 Current address: Department of Medical Physics, Memorial SloanKettering Cancer Center, New York, NY 10021. 2 Current address: Department of Radiology, Division of Neuroradiology, University of Texas Southwestern Medical Center at Dallas, Dallas, TX 75390. Available online on ScienceDirect (www.sciencedirect.com.) 1053-8119/$ - see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2004.03.024

duration, inter-stimulus interval duration, and incompletely understood interrelationships among neural activity, cerebral metabolic rate of oxygen consumption (CMRO2), cerebral blood flow (CBF), and cerebral blood volume (CBV) (Buxton et al., 1998). If the stimulus trials are repeated identically and separated sufficiently in time, the experimenter can assume that the hemodynamic response will be the same each time in activated cortical regions (Aguirre et al., 1998; Miezin et al., 2000). However, if a stimulus presentation parameter changes, the experimenter cannot make this assumption, but must know the relationship between the stimulus parameter and the hemodynamic response to accurately interpret the functional data. Early experiments (Boynton et al., 1996) pointed to a linear relationship between stimulus parameters, such as duration, and the hemodynamic response. Others (Cohen, 1997; Friston et al., 1994) have shown that a linear systems analysis can be used to model the BOLD response. However, more recent papers and abstracts (Table 1) have shown that nonlinearity does indeed exist between the stimulus parameters (duration, rate, and amplitude) and the resulting hemodynamic response for various types of stimuli. The origin of the nonlinearity of the hemodynamic response remains a mystery. Since a cascade of physiological processes occurs after stimulus presentation, it is likely that multiple factors between the stimulus and the measured BOLD response contribute to nonlinearity. Nonlinearities in the responses of neuronal activity, oxygen extraction, blood flow, and blood volume are possible, for example. Rees et al. (1997) found a nonlinear BOLD response and a linear CBF response in the auditory cortex as a function of presentation rate, suggesting that a factor other than CBF is responsible for the nonlinearity. The authors suggested that the nonlinearity might be caused by changes in the oxygen extraction fraction (OEF). However, Miller et al. (2001) reported a strongly nonlinear response as a function of stimulus duration for both CBF and BOLD in the visual cortex. Logethetis et al. (2001) reported a linear relationship between local field potentials (LFP) of neurons and the BOLD response and stated that both LFPs and the BOLD response were nonlinear functions of visual stimulus contrast. This result suggests that any nonlinearity of the BOLD

1118

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

Table 1 Papers and abstracts showing nonlinearity of the BOLD response Stimulus, Cortex

Visual

Motor

Auditory

Stimulus parameter Duration

Rate

Amplitude

Pfeuffer et al., 2003 Birn et al., 2001 Miller et al., 2001 Liu and Gao, 2000 Vazquez and Noll, 1998 Birn et al., 2001 Miller et al., 2001 Glover, 1999 Glover, 1999

Singh et al., 2003 Chawla et al., 1999 Kwong et al., 1992

Vazquez and Noll, 1998

Ja¨ncke et al., 1998

Peck et al., 2001

Robson et al., 1998

Somatosensory

Nangini et al., 2002

(Language Task: Multiple Cortices) (test with additivity)

Friston et al., 1998 Rees et al., 1997 Binder et al., 1994 Harrington and Downs, 2001 Takanashi et al., 2000 Gopinath et al., 2001 Mechelli et al., 2000 (test with scaling)

(test with scaling)

response may be due to nonlinearities in the neural response. Although one study (Birn et al., 2001) did not find any correlation between response latency and degree of nonlinearity, another study (Pfeuffer et al., 2003) did find such a correlation, suggesting that vessel size may be a contributing factor to nonlinearity. The exact origin of the nonlinear effect is still not clear. Perhaps, the nonlinearity of the BOLD response is caused by a combination of factors, such as neural adaptation, blood flow, and oxygen extraction. Few of the studies represented in Table 1 describe the ranges of linearity or nonlinearity for the BOLD response. Robson et al. (1998) reported that the BOLD response in the primary auditory cortex is nonlinear for stimuli less than 6 s and linear above this duration. Two studies looking at the primary visual cortex reported that the BOLD response is nonlinear for stimuli less than 4 s (Vazquez and Noll, 1998) or nonlinear for stimuli less than 3 s (Liu and Gao, 2000), and linear above these durations. No one has reported such ranges for the primary motor cortex. Knowledge of these ranges for different cortices would be very helpful in designing experiments where stimulus duration is likely to change. Two studies (Birn et al., 2001; Pfeuffer et al., 2003) have looked at the variation of nonlinearity on a spatial basis. Individual voxels were assigned indices of nonlinearity based on the response profiles for several stimulus or task durations. In a motor task experiment, Birn et al. (2001) found that the measure of nonlinearity for voxels differed between the SMA and primary motor

cortex. This study also found voxelwise variation of the nonlinearity index in the primary visual cortex, but no clear difference between the visual and motor cortices. Pfeuffer et al. (2003) looked at the spatial dependence of the nonlinearity only within the visual cortex. Therefore, further research is necessary to compare the nonlinearity across different cortices of the brain. A linear system is defined by the principle of superposition (Oppenheim et al., 1999). If the responses of a system to inputs x1 and x2 are T(x1) and T(x2), then a linear system requires T(x1 + x2) is equal to T(x1) + T(x2). This is the property of additivity. In addition, if the response of a system to input a
x, where a is a scalar constant, is T(ax) then a linear system requires T(ax) is equal to aT(x). This is the property of scaling. It is important to note that increasing the duration of the stimulus will not cause a scaled increase in the resultant hemodynamic response in a linear system as it would for a change in stimulus intensity. The effect of changing duration can be predicted with a test of additivity, that is, summing the responses to shorter stimuli, the durations of which add up to the duration of the longer stimulus. For the purposes of this experiment, only the parameter of duration was tested in examining the linearity of the BOLD response. Another way to test the linearity of a system is by trying to estimate the impulse response function of the response due to several different inputs (Liu and Gao, 2000). The impulse response function in a discrete time system can be defined as the response of a system to an input stimulus of unit time, or 1 s in the case of typical fMRI data. The response to a stimulus of greater duration can be found by convolving the stimulus function, x(t), with the impulse response function, i(t), rðtÞ ¼ xðtÞ  iðtÞ ¼

l X

xðtÞiðt  kÞ

ð1Þ

k¼0

In a linear system, the impulse response function would be the same for any stimulus duration. In a nonlinear system, the impulse response function would not be the same for different stimulus durations. Knowing how the shape of the impulse response function changes with changes in stimulus duration can help to visualize the point at which the system changes from linear to nonlinear. The main purpose of this study was to examine the ranges of linearity and nonlinearity for the relationship between stimulus duration and hemodynamic response for three primary cortices (auditory, motor, and visual). A single group of subjects was used to examine any differences that might exist in the nonlinearity between cortices. The results of this analysis will be useful for designing fMRI experiments where the stimulus or task duration is not constant. In addition, the results of two different methods for examining the nonlinearity of the BOLD response with respect to duration were compared. Finally, the nonlinearity of the hemodynamic response was modeled using a neural adaptation function.

Materials and methods Subjects Five right-handed male volunteers (19.9 F 1.9 years old) were scanned. The subjects gave informed consent according to local IRB protocol. For a fixed-effect analysis, the use of five subjects allowed an inference of 55% of the population with a = 0.05 and b = 0.05 (Friston et al., 1999). The subjects were age- and gender-

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

1119

Fig. 1. Averaged HRFs for the primary auditory, motor, and visual cortices in one subject. SEM error bars are barely visible. Colored horizontal bars indicate the duration of the various stimuli.

matched to avoid variations in the BOLD response due to differences caused by age (D’Esposito et al., 1999) and gender (Levin et al., 1998). Each subject reported having no hearing impairment, corrected vision, or known history of hypertension.

white cycle subtended about 1.1j. The control condition was a gray screen of equal luminance to the checkerboard. This type of visual stimulus was similar to that used by many researchers (Birn et al., 2001; Boynton et al., 1996; Vazquez and Noll, 1998).

Stimuli/task

Experiment design

Two types of stimuli (auditory and visual) and one type of task (manual motor) were used in this experiment. The auditory stimuli consisted of a train of 440 Hz tones, each lasting 125 ms followed by a rest period of 125 ms, producing four tones per second. The resting state was silence. The motor task, run at the same time as the auditory stimulus, consisted of the subject bringing four fingers of the right hand together with the thumb of the same hand four times a second in unison with the pulsating tones. The resting condition was no motion. The frequency of the motor task was kept near the maximum possible rate, as that is known to elicit strong motor cortex activation (Ja¨ncke et al., 1998). The combined nature of the auditory stimulus and motor task is similar to that used by Glover (1999). The visual stimulus consisted of a contrast-reversing black and white checkerboard alternating at a frequency of 7.5 Hz, as it is known to elicit strong visual cortex activation (Fox and Raichle, 1984; Kwong et al., 1992; Singh et al., 2003). The 24  16 checkerboard filled a visual field of 11j  7.3j and each black and

Two separate scanning sessions were run on each subject. The first consisted of the auditory stimuli and motor tasks executed simultaneously. The visual stimulus experiment was run on the second day. The experimental design for each scanning session was identical. The first two functional runs consisted of 10 stimulus cycles. Each cycle consisted of 12 s of stimulus or task followed by 12 s of rest in a basic block design. Activated voxels would be chosen from these runs. Eight functional runs followed, each consisting of 10 stimulus cycles. Each cycle had one of five different duration periods (1, 2, 4, 8, or 16 s) chosen in a pseudorandom fashion followed by 20 s of the resting condition. The BOLD responses for each of the different stimulus duration periods were obtained from these eight runs. Data acquisition All data were acquired on a GE 3T LX scanner (General Electric Medical Systems, Milwaukee, WI) using a dome-shaped

Fig. 3. Idealized impulse response functions for stimuli of five different duration periods for the primary auditory, motor, and visual cortices in one subject.

1120

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

quadrature RF head coil (MRI Devices, Waukesha, WI). A T1weighted spoiled GRASS sequence was used to obtain highresolution anatomical images, and a time-of-flight MRA was used to map blood vessels. Functional data were acquired using a 1shot spiral sequence (TR/TE/FA = 1000 ms/18 ms/60j) (Noll et al., 1995). Twenty-three-millimeter-thick coronal slices were obtained with a FOV of 200 mm and a matrix size of 64  64, yielding an in-plane resolution of 3.125  3.125 mm2. The slices were carefully prescribed to cover the cortex of interest by consulting a brain atlas (Mai et al., 1997). The entire brain was not scanned to improve the temporal resolution of the imaging sequence. Data analysis All data were analyzed using AFNI (Cox, 1996) and MATLAB (The Mathworks, Natick, MA). Three dimensional volume registration was applied to all functional image time series (Cox and Jesmanowicz, 1997). Linear trends were removed before separate runs were concatenated. The first two runs, consisting of periodic stimuli or motor performance (12 s on/12 s off) were concatenated into what shall be called Run 1, which was used to determine the activated voxels in the data set. The next eight runs, called Runs 2a – h, were used to characterize the shape of the hemodynamic response for stimuli of different duration. A reference vector was constructed, consisting of the positive lobe of a sine curve (with a period equal to 24 s) for the first 12 images acquired during the stimulus or motor performance and zero for the 12 images acquired during rest. This cycle was repeated 20 times to match the length of Run 1, and 24 timeshifted copies were made. An initial cross correlation analysis (Bandettini et al., 1993) was done between these reference vectors and each of the voxel time series in the data set of Run 1. The reference vector that gave the highest positive correlation value was assumed to have the correct delay for the hemodynamic response. The time series of the voxel with the highest positive correlation, using the appropriately time-shifted reference vector,

was then averaged across all 20 cycles and smoothed with a 3neighborhood order-statistics filter (Bovik et al., 1983), yn ¼ 0:15 minða; b; cÞ þ 0:70medianða; b; cÞ þ 0:15 max ða; b; cÞ;

ð2Þ

to create the time-averaged response. This time-averaged response was then convolved with a stimulus vector representing a value of 1 for the onset of each of the 20 cycles of stimulus to create the time-averaged response vector. It was assumed that this timeaveraged response vector would contain the appropriate delay and width of the hemodynamic response and thus provide a more accurate model for the HRF than an arbitrarily modeled waveform. To assure that the model HRF would contain the appropriate temporal characteristics for each particular region, separate timeaveraged response vectors were made for each of the cortices (auditory, motor, and visual). The three empirically derived reference vectors were then used in a second cross correlation analysis to mask the voxels used for analysis of Run 2 data. Voxels with a correlation coefficient (r) greater than 0.5 ( P < 6.6  1027 Bonferroni corrected) were defined to be active. For each experiment, activated voxels were identified as being in the primary auditory, motor, or visual cortex with the aid of an anatomical brain atlas (Mai et al., 1997). The eight runs containing randomized stimuli of the five different duration periods (Run 2a – h) were sorted to create five separate time course data sets, one for each duration. The time series data corresponding to the trials for the same stimulus duration were concatenated to produce separate fMRI time series. Thus, five data sets were obtained, each containing 16 trials for one specific duration period. Temporal smoothing was applied to these time courses with a 3neighborhood order-statistics filter, as given in Eq. (2). The hemodynamic responses were then averaged across all 16 trials and all activated voxels found in the appropriate cortex in the block design experiment of Run 1. In the auditory cortex, subjects activated volumes ranging from 0.94 (32 voxels) to 2.23 cc (76

Fig. 2. The peak amplitude from the predicted responses are shown alongside the 95% confidence intervals of the observed responses for the 16-s responses in the primary auditory, motor, and visual cortices. The abscissa refers to the stimulus duration that produced the response used to predict the response for the 16-s stimulus. Confidence intervals and standard deviation error bars were determined from the responses of the five subjects.

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127 Table 2 Impulse response function parameters

1121

Results

Stimulus duration (s) Amplitude (%) Area (% s) TTP (s)

FWHM (s)

Averaged HRFs

Auditory cortex 1 2 4 8 16

2.1 1.6 1.2 0.8 0.5

F F F F F

0.3 0.1 0.1 0.1 0.1

7.1 6.2 4.6 3.2 2.9

F F F F F

1.6 0.7 0.9 1.0 1.0

4.1 4.1 3.7 3.1 3.0

F F F F F

0.2 0.3 0.2 0.5 0.8

3.2 3.6 3.5 3.8 4.9

F F F F F

0.4 0.5 0.5 0.9 0.7

Motor cortex 1 2 4 8 16

2.0 1.6 1.2 0.7 0.7

F F F F F

0.8 0.4 0.3 0.2 0.2

8.8 8.0 6.1 3.7 3.9

F F F F F

4.3 2.9 1.6 1.1 0.7

5.1 5.0 4.9 4.4 4.7

F F F F F

0.6 0.6 0.3 0.6 0.3

4.1 4.5 4.7 4.6 5.4

F F F F F

0.7 0.6 0.6 0.5 1.2

Visual cortex 1 2 4 8 16

2.2 1.6 1.2 1.1 1.2

F F F F F

0.4 0.3 0.2 0.1 0.2

10.5 7.7 6.6 5.8 5.5

F F F F F

3.1 1.9 1.4 0.8 0.6

4.3 4.4 4.8 5.0 5.5

F F F F F

0.5 0.6 0.7 0.5 0.6

4.3 4.5 5.1 5.0 4.1

F F F F F

0.6 0.7 0.7 0.8 0.4

Results for the hemodynamic responses in the primary auditory, motor, and visual cortices activated for the five different duration periods in one subject are shown in Fig. 1. The amplitudes for the auditory and motor responses reach a plateau after roughly 4 s of stimulus or task presentation. The amplitude for the visual response appears to reach a plateau at a point between 8 and 16 s of stimulus presentation. The response amplitudes appear to match across all three cortices for the shorter stimuli, but vary for the longer stimuli. An analysis of variance using responses across subjects was performed to determine the significance of the difference in amplitude across cortices. A test of equality was done for the amplitudes of the responses across cortices. The F-statistics for the 8- and 16-s responses were 16.77 and 22.88, respectively, both of which are significant for P < 0.001. Thus, at least one cortex produced a population of amplitudes that was different from the other cortices. No significant differences were found for the amplitudes of the 1-, 2-, and 4-s responses. Additivity method

voxels). In the motor cortex, subjects activated volumes ranging from 1.41 (36 voxels) to 3.43 cc (117 voxels). In the visual cortex, subjects activated volumes ranging from 1.64 (56 voxels) to 4.66 cc (159 voxels). Averaging a large number of responses yields a high signal-to-noise ratio (SNR) for the response, assuming intraregional homogeneity. To find the percent signal change of the hemodynamic response, a baseline for each voxel time series had to be determined. The first value of each trial, corresponding to the start of each stimulus, was averaged across all 16 trials to produce the baseline. This was done separately for each time course corresponding to a different duration period. It was assumed that the rest period of 20 s, following each stimulus event, would allow the hemodynamic response enough time to return to baseline (Dale and Buckner, 1997).

The additivity test was first used to examine the linearity of the BOLD response. The HRFs for the 1-, 2-, 4-, and 8-s stimuli were used to predict the responses to the 2-, 4-, 8-, and 16-s stimuli. For example, the 1-s response was shifted by 1 s and added to the 1-s response to predict the HRF to the 2-s stimulus. Over-predictions were apparent in most cases, most especially seen in peak amplitude differences. To assess the significance of the differences in the amplitude between the observed and predicted responses, 95% confidence intervals for the amplitudes of the observed responses were compared with the predicted amplitude responses. In Fig. 2, the predictions for the 16-s response amplitude based on the responses for the 1-, 2-, 4-, and 8-s responses are examined. For most cases, the prediction of the 16-s response amplitude falls

Fig. 4. The amplitude of the impulse response function decreases with increasing stimulus duration across subjects. The 95% confidence intervals for the amplitude of the impulse function to the 16-s stimulus are shown. From these curves, one can estimate the ranges of linearity and nonlinearity for the BOLD response in the primary auditory, motor, and visual cortices. Confidence intervals and standard deviation error bars were determined from the impulse response functions of the five subjects.

1122

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

Fig. 5. Multiplying the normal stimulus function, s(t), with the neural adaptation function, n(t), yields a modified stimulus function, which can then be convolved with the impulse response function, i(t) to yield the BOLD response, r(t).

above the confidence intervals for the observed 16-s response amplitude. From the intersection of the predicted response curve with the 95% confidence interval of the observed 16-s response amplitude, the BOLD response can be seen to change from nonlinear to linear at a certain threshold. This threshold is roughly estimated to be 7 s for the auditory cortex, 8 s for the motor cortex, and indeterminate for the visual cortex. Impulse response function method A second test of the linearity was performed by finding the impulse response function for each response, using time (s) for the abscissa and percent change (%) for the ordinate. An algorithm was developed to fit the impulse response function to a three-parameter gamma – variate function, iðtÞ ¼ c1 ðt=c2 Þc3 expðt=c2 Þ;

ð2Þ

and then convolve it with a rectangular stimulus function (of height equal to 1 and duration equal to the stimulus duration) to produce a fitted response function. Ranges of values were tried for the

parameters c1, c2, and c3. The sum of squares between the observed HRF and the fitted response function was minimized by finding the best values for c1 (which affects the amplitude of the response function), c2 (which affects the full width at half maximum (FHWM) and the time to peak (TTP) of the response function), and c3 (which affects the amplitude and TTP of the response function). The algorithm normalized the fitted response function and searched for a single constant at a time (c2, c3, and finally c1) to reduce search time. Although the fitted HRF fails to predict the undershoot or the double peak, it does fit the post-stimulus overshoot remarkably well. Such a model assumes a linearity of the BOLD response. The impulse response function parameters averaged across subjects are shown in Table 2. If the BOLD response was linear, the values for each parameter would be consistent, but an examination of the values for amplitude and area clearly show significant variation. The amplitude for the impulse response functions for the 1-s stimulus is similar across all three cortices. However, for the 16s stimulus, the amplitude for the impulse response function in the visual cortex is larger than that of the auditory or motor cortex.

Fig. 6. (A) The impulse response functions are seen alongside (B) the neural adaptation functions for the primary auditory, motor, and visual cortices in all five subjects.

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

1123

Table 3 Parameters for the impulse response function and the neural adaptation function

Auditory cortex Motor cortex Visual cortex

Amplitude

Area

TTP

FWHM

a

b

0.7 F 0.2 0.6 F 0.1 0.7 F 0.1

2.7 F 0.9 3.4 F 0.7 5.4 F 0.5

4.7 F 0.3 5.8 F 0.6 5.3 F 0.7

3.8 F 0.5 5.7 F 1.3 7.4 F 0.3

2.4 F 1.5 2.2 F 0.7 1.5 F 0.6

0.7 F 0.4 0.7 F 0.4 0.6 F 0.4

A graphical representation of the variation in the impulse response function across stimulus duration periods is shown in Fig. 3. As the stimulus duration increases, the amplitude of the impulse response function decreases. The impulse response function amplitude is plotted as a function of stimulus duration for all three cortices in Fig. 4, along with 95% confidence intervals for the amplitude of the impulse function for the 16-s stimulus (which is assumed to be in the range of linearity). The transition from nonlinearity to linearity of the BOLD response can be defined as the point (stimulus duration) for which the curve representing the amplitude of the impulse response function as a function of stimulus duration crosses into the confidence interval region. The threshold can be estimated to be about 10 s for the auditory cortex, 7 s for the motor cortex, and 3 s for the visual cortex. Besides amplitude, other factors of the impulse response function were compared, including area, FWHM, and TTP. Significant trends were seen for the area as a function of stimulus duration as would be expected from the significant trends seen for amplitude. However, significant trends were not seen for the FWHM or the TTP of the impulse response function as a function of stimulus duration. Neural adaptation A pair of studies (Miller et al., 2001; Pfeuffer et al., 2003) have tried to explain the nonlinearity of the hemodynamic response as a result of neural adaptation. Neural adaptation can be modeled by a simple function, nðtÞ ¼ 1 þ aexpðt=bÞ

ð4Þ

where a represents the initial offset for the level of neuronal firing and b represents the decay constant. The decrease in amplitude of the impulse response function with increased stimulus or task duration (Fig. 3) can be interpreted as being due to neural adaptation, which refers to the decay of neural spike rates for a sustained stimulus (Rieke et al., 1997). In reality, the nonlinearity of the hemodynamic response may be a result of several factors, including neural adaptation as well as changes in cerebral blood flow and oxygen extraction fraction. The product of the stimulus function and the neural adaptation function would then be convolved with the impulse response function to obtain the resulting hemodynamic response (Fig. 5). In this way, a nonlinear model would be used to explain the nonlinearity of the hemodynamic response. For each subject, our data was analyzed to solve for a single impulse response function and a neural adaptation function for all durations. For this algorithm, a pair of nested loops was used to solve for the best values of a and b in Eq. (4) using a sum of squares minimization method, comparing the observed responses with the predicted responses. At the same time, the values for c1, c2, and c3 for the impulse response function were found using the algorithm described earlier, but applied simultaneously to all five responses. Since the 16-s responses were observed with more time

points than the 1-s responses, a weighting factor was used to make sure that each response was treated equally. The resulting impulse response functions and neural adaptation functions for each subject and each primary cortex are shown in Fig. 6. As can be seen, within a cortex, the neural adaptation functions of different individuals are fairly similar with one another with a few exceptions. Table 3 shows the average impulse response function parameters and average values of the intensity overshoot a and the exponential time constant b from Eq. (4). The impulse response function amplitude appears quite consistent across the three cortices for this model. For the motor cortex, Miller et al. (2001) reported values of a = 0.25 and b = 0.25, while we found a = 2.2 F 0.7 and b = 0.7 F 0.4. The auditory cortex showed similar values with a = 2.4 F 1.5 and b = 0.7 F 0.4. For the visual cortex, Miller et al. (2001) reported a = 3 and b = 0.5, while we found a = 1.5 F 0.6 and b = 0.6 F 0.4. The performance of the linear model and the nonlinear model (with the neural adaptation function) were compared to assess any significant improvement afforded by the latter. To do this, a single impulse response function was estimated based on the averaged BOLD responses to the five different stimulus/task duration periods for the linear model just as was done for the nonlinear model as described above. Then, the sum of squares between the fitted responses and the observed responses were found for both the linear and nonlinear model. It was found that in all cases, the sum of squares for the nonlinear model was lower than that of the linear model (Table 4). Furthermore, F-tests

Table 4 Comparison of linear and nonlinear models Subject

Sum of squares Linear model

F(2,126)

P value

Nonlinear model

Auditory cortex Sub01 98.58 Sub02 38.11 Sub03 35.01 Sub04 35.92 Sub05 26.94

45.54 18.40 22.50 21.31 17.00

73.24 67.49 35.03 43.19 36.84

7.89 1.19 8.00 5.19 2.53

    

1022 1020 1013 1015 1013

Motor cortex Sub01 133.09 Sub02 21.50 Sub03 64.14 Sub04 41.45 Sub05 51.18

54.53 11.43 41.08 12.33 19.14

90.76 55.50 35.36 148.79 105.46

3.87 5.18 6.47 6.71 1.23

    

1025 1018 1013 1034 1027

Visual cortex Sub01 23.97 Sub02 133.76 Sub03 34.33 Sub04 134.97 Sub05 52.05

15.67 19.67 19.03 28.32 17.58

33.37 365.41 50.65 237.25 123.53

2.34 3.56 7.20 1.89 2.00

    

1012 1053 1017 1043 1030

1124

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

were performed to assess the significance in the reduction of the sum of squares after adding two neural adaptation parameters to the model, where: FðDF1  DF2; DF2Þ ¼

ðSS1  SS2Þ=SS2 : ðDF1  DF2Þ=DF2

ð5Þ

The degrees of freedom for each model were calculated as the number of time points for all five responses (131) minus the number of parameters in the model (3 for the linear model, 5 for the nonlinear model). In all cases, the nonlinear model produced highly significant reductions in the unexplained variance when compared to the linear model (Table 4). One interesting result of using a neural adaptation function is that the shape of the fitted hemodynamic response is able to model a peak followed by a plateau for the 16-s response (Fig. 7B). The convolution of a rectangular stimulus function and an impulse response function modeled as a gamma-variate results only in a simple response peak (Fig. 7A). It thus appears that the concept of neural adaptation might explain the non-intuitively shaped peak sometimes seen for BOLD responses due to long stimuli. It should be noted that this biphasic response was only seen for the 16-s response in three out of five subjects for the auditory cortex, two

out of five subjects for the motor cortex, and in no subjects for the visual cortex.

Discussion This study examined the nonlinearity of the BOLD response with respect to duration in three primary cortices. Ranges of linearity and nonlinearity were assessed for each cortex. A single group of subjects was used allowing the comparison of these ranges across primary cortices. In addition, a neural adaptation model was used to explain the observed nonlinearity. The two separate methods of testing the linearity of the hemodynamic response produced somewhat different results (Table 5). In the primary auditory cortex, the BOLD response was seen to be nonlinear for stimuli of duration less than about 7 s for the additivity method and nonlinear for stimuli less than 10 s for the impulse response function method. The second result is somewhat different than that found by Robson et al. (1998), who found that the nonlinear range in the auditory cortex occurs for stimuli less than 6 s. The difference in nonlinearity analysis methods could account for this discrepancy. In the primary motor cortex, the BOLD response was seen to be nonlinear for stimuli of duration

Fig. 7. The observed 16-s responses (circles) and the fitted responses (solid lines) in the primary auditory, motor, and visual cortices in one subject. (A) The fitted responses were obtained by convolving a rectangular stimulus function with the idealized impulse response function, represented by a gamma-variate. (B) The fitted responses were obtained by modeling a neural adaptation function along with the impulse response function. Shown are the correlation coefficients (r) between the observed and fitted values for each case.

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127 Table 5 Ranges of linearity and nonlinearity with respect to duration Auditory cortex (s) Additivity method Nonlinear range Linear range

<7 >7

Impulse response function method Nonlinear range <10 Linear range >10

Motor cortex (s)

Visual cortex (s)

<8 >8

V8 not determined

<7 >7

<3 >3

less than 7 or 8 s and linear for stimuli above this threshold for both tests. No report is made in the literature for ranges of nonlinearity in the motor cortex. In the primary visual cortex, the BOLD response was seen to be nonlinear for stimuli of duration less than at least 8 s for the additivity method and nonlinear for stimuli less than 3 s for the impulse function method. The second result agrees precisely with that found by Liu and Gao (2000). Liu and Gao (2000) also used an impulse response function method to assess the linearity of the system. Vazquez and Noll (1998) found that the nonlinear range in the visual cortex occurs for stimuli less than 4 s. The method of additivity is simple to execute, but could include errors resulting from a poor estimation of the complete hemodynamic responses. After averaging the responses from many trials and many voxels, the final HRF shape is likely to contain deviations (particularly at the tail end of the response), which will confound the accurate comparison of predicted and observed responses. The undershoot of the hemodynamic response is a small change from the baseline (approximately 0.5 to 1%), and thus is difficult to separate from the noise of the time series, which has a standard deviation about 0.5% to 1% of the baseline as seen in our data. The second method, which determines an idealized estimate of the impulse response function, avoids the problem of the undershoot by ignoring it completely. Thus, this test of linearity concentrates solely on the post-stimulus overshoot, which is the main part of the BOLD response. For this reason, the results of the impulse response function method seen in Table 5 are different from those of the additivity method. Since the detection of activation commonly relies on the BOLD overshoot, it makes sense to restrict the linearity analysis to this part of the hemodynamic response. Neural adaptation has been shown as a convenient way to account for the nonlinearity of the BOLD response for variation in stimulus duration. However, neural adaptation is not necessarily the sole cause of the nonlinearity. Other physiological processes may also contribute to this phenomenon. Miller et al. (2001) indicated that the BOLD response saturates at high levels of cerebral blood flow. This means that using BOLD responses during low levels of CBF would result in over-predictions for BOLD responses at high levels of CBF. Buxton and Frank (1997) revealed that the intensity of the fMRI signal is strongly coupled to the oxygen extraction fraction. If the OEF changes during a sustained stimulus, the resultant BOLD signal will be affected. Several potential sources of error exist for this study. It was assumed that all of the voxels in a particular region would produce similar BOLD responses, as reported by others (Miezin et al., 2000). However, it is possible that this is not true, as some researchers (Birn et al., 2001; Pfeuffer et al., 2003) have reported

1125

a variation of nonlinearity on a voxelwise basis. However, such studies suffer from a limited number of trials, thus yielding a low SNR for the averaged hemodynamic response. In addition, it has been reported that the temporal delay of the BOLD response varies as a function of blood vessel diameter (Krings et al., 1999; Lee et al., 1995). Thus, if differently sized blood vessels contributed to the hemodynamic responses in different activated voxels, the shape of the averaged response would be affected. Methodological sources of error could also exist. It was assumed that a rest period of 20 s would allow the hemodynamic signal to return to baseline. If the signal does not fully recover by the time of the new stimulus, the baseline estimate will be affected, leading to incorrect values for percent signal change. It should be noted that each run contained trials of each stimulus duration period arranged in a pseudo-random order. So, if the longer duration stimuli produced a more pronounced post-stimulus undershoot that did not recover within the 20-s rest period, responses of all other duration periods would have been affected. Longer rest periods could be used to ensure a proper return to baseline, but this would either increase the length of the experiment or reduce the number of trials used for each duration period. In addition, there could be low frequency fluctuations present in the data from scanner instability that would confound the estimate of the baseline.

Conclusion The research on the nonlinearity of the hemodynamic response presented in this study has yielded ranges of stimuli or task duration for which the BOLD response is nonlinear or linear with respect to duration. The ranges of nonlinearity have been found to vary between cortices. This variation may be caused by regional differences in neuronal output or cerebral vasculature. Researchers desiring to stimulate cortical areas of the brain for varying amounts of time may want to take into account the problems associated with the nonlinearity of the hemodynamic response. Specifically, if subjects in the magnet are exposed to stimuli or engaged in tasks for both short- and long-duration periods, a linear model may not accurately predict the shape of the hemodynamic response for all such duration periods. To avoid any potential problems associated with the nonlinearity, researchers may want to restrict stimuli or task periods to the range of linearity for the cortex to be studied. The neural adaptation functions could offer a way to correct for the nonlinearity of the hemodynamic response with respect to duration. To locate activated voxels in a typical fMRI experiment, researchers could create a reference wave that incorporates the nonlinear model by convolving an appropriate impulse response function with a stimulus duration function that has been multiplied by an appropriate neural adaptation function. The parameters for the neural adaptation function obtained in this experiment could be used for such analysis. However, variability in scanners, subjects, and stimulus intensity are likely to cause variability in these parameters. Additional research is needed to examine how the detection of activation would change with the use of neural adaptation functions. In addition, multimodal experiments could be done to compare the level of EEG activity with the neural adaptation functions derived in these fMRI experiments. Such work would verify if a

1126

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127

major source of the nonlinearity is neural adaptation, or a decay of neural spike rates for a sustained stimulus (Rieke et al., 1997), and provide more accurate data for further modeling of the hemodynamic responses observed in fMRI. One study of visual evoked potentials (VEP) and BOLD responses (Janz et al., 2001) suggested that neural adaptation alone cannot explain the nonlinear behavior of hemodynamic response. The presence of nonlinearity in the hemodynamic response reminds us that more research must be done to understand the exact nature of the relationship between the neural activity and the BOLD response observed in fMRI.

Acknowledgments This research was supported in part by NIH grant P50DC03888, the Center of Excellence Grant #F2182C, and the VA RR&D Brain Rehabilitation Center, Gainesville, FL.

References Aguirre, G.K., Zarahn, E., D’Esposito, M., 1998. The variability of human, BOLD hemodynamic responses. NeuroImage 8, 360 – 369. Bandettini, P.A., Jesmanowicz, A., Wong, E.C., Hyde, J.S., 1993. Processing strategies for time-course data sets in functional MRI of the human brain. Magn. Reson. Med. 30, 161 – 173. Birn, R.M., Saad, Z.S., Bandettini, P.A., 2001. Spatial heterogeneity of the nonlinear dynamics in the FMRI BOLD response. NeuroImage 14, 817 – 826. Binder, J.R., Rao, S.M., Hammeke, T.A., Frost, J.A., Bandettini, P.A., Hyde, J.S., 1994. Effects of stimulus rate on signal response during functional magnetic resonance imaging of auditory cortex. Cognit. Brain Res. 2, 31 – 38. Bovik, A.C., Huang, T.S., Munson, D.C., 1983. A generalization of median filtering using linear combinations of order statistics. IEEE Trans. Acoust., Speech, Signal Process. 31, 1342 – 1349. 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. Buxton, R.B., Frank, L.R., 1997. A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation. J. Cereb. Blood Flow Metab. 17 (1), 64 – 72. Buxton, R.B., Wong, E.C., Frank, L.R., 1998. Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magn. Reson. Med. 39, 855 – 864. Chawla, D., Buechel, C., Edwards, R., Howseman, A., Josephs, O., Ashburner, J., Friston, K.J., 1999. Speed-dependent responses in V5: a replication study. NeuroImage 9, 508 – 515. Cohen, M.S., 1997. Parametric analysis of fMRI data using linear systems methods. NeuroImage 6, 93 – 103. Cox, R.W., 1996. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput. Biomed. Res. 29, 162 – 173. Cox, R.W., Jesmanowicz, A., 1997. Real-time 3D image registration for functional MRI. Magn. Reson. Med. 42, 1014 – 1018. Dale, A.M., Buckner, R.L., 1997. Selective averaging of rapidly presented individual trials using fMRI. Hum. Brain Mapp. 5, 329 – 340. D’Esposito, M., Zarahn, E., Aguirre, G.K., Rypma, B., 1999. The effect of normal aging on the coupling of neural activity to the BOLD hemodynamic response. NeuroImage 10 (1), 6 – 14. Fox, P.T., Raichle, M.E., 1984. Stimulus rate dependence of regional cerebral blood flow in human striate cortex, demonstrated by positron emission tomography. J. Neurophysiol. 51, 1109 – 1120. Friston, K.J., Jezzard, P., Turner, R., 1994. Analysis of functional MRI time-series. Hum. Brain Mapp. 1, 153 – 171.

Friston, K.J., Josephs, O., Rees, G., Turner, R., 1998. Nonlinear eventrelated responses in fMRI. Magn. Reson. Med. 39, 41 – 52. Friston, K.J., Holmes, A.P., Price, C.J., Buchel, C., Worsley, K.J., 1999. Multisubject fMRI studies and conjunction analysis. NeuroImage 10, 385 – 396. Glover, G.H., 1999. Deconvolution of impulse response in event-related BOLD fMRI. NeuroImage 9, 416 – 429. Gopinath, K.S., Briggs, R.W., Himes, N., 2001. Examination of the linearity of BOLD FMRI responses in a higher level cognitive system. Proceedings of 9th ISMRM, Glasgow, p. 1189. Harrington, G.S., Downs III, J.H., 2001. FMRI mapping of the somatosensory cortex with vibratory stimuli: is there a dependency on stimulus frequency? Brain Res. 897, 188 – 192. Hennig, J., Speck, O., Koch, M.A., Weiller, C., 2003. Functional magnetic resonance imaging: a review of methodological aspects and clinical applications. JMRI 18, 1 – 15. Ja¨ncke, L., Specht, K., Mirzazade, S., Loose, R., Himmelbach, M., Lutz, K., Shah, N.J., 1998. A parametric analysis of the ‘rate effect’ in the sensorimotor cortex: a functional magnetic resonance imaging analysis in human subjects. Neurosci. Lett. 252, 37 – 40. Janz, C., Heinrich, S.P., Kornmayer, J., 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. Magn. Reson. Med. 46, 482 – 486. Krings, T., Erberich, S.G., Roessler, F., Reul, J., Thron, A., 1999. MR blood oxygenation level-dependent signal differences in parenchymal and large draining vessels: implications for functional MR imaging. Am. J. Neuroradiol. 20, 1907 – 1914. 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. Lee, A.T., Glover, G.H., Meyer, C.H., 1995. Discrimination of large venous vessels in time-course spiral blood – oxygen-level-dependent magneticresonance functional neuroimaging. Magn. Reson. Med. 33, 745 – 754. Levin, J.M., Ross, M.H., Mendelson, J.H., Mello, N.K., Cohen, B.M., Renshaw, P.F., 1998. Sex differences in blood-oxygenation-level-dependent functional MRI with primary visual stimulation. Am. J. Psychiatry 155, 434 – 436. Liu, H.L., Gao, J.H., 2000. An investigation of the impulse functions for the nonlinear BOLD response in functional MRI. Magn. Reson. Med. 18, 931 – 938. Logethetis, N.K., Pauls, J., Augath, M., Trinath, T., Oeltermann, A., 2001. Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 150 – 157. Mai, J.K., Assheuer, J., Paxinos, G., 1997. Atlas of the Human Brain. Academic Press, San Francisco, CA. Mechelli, A., Friston, K.J., Price, C.J., 2000. The effects of presentation rate during word and pseudoword reading: a comparison of PET and fMRI. J. Cogn. Neurosci. 12 (Suppl. 2), 145 – 156. Miezin, F.M., Macotta, L., 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. Miller, K.L., Luh, W.M., Liu, T.T., Martinez, A., Obata, T., Wong, E.C., Frank, L.R., Buxton, R.B., 2001. Nonlinear temporal dynamics of the cerebral blood flow response. Hum. Brain Mapp. 13, 1 – 12. Nangini, C., Staines, W.R., McIlroy, W.E., Graham, S.J., 2002. Non-linear event related fMRI BOLD responses in human somatosensory cortex. Proceedings of 10th ISMRM, Honolulu, p. 1385. Noll, D.C., Cohen, J.D., Meyer, C.H., Schneider, W., 1995. Spiral k-space MR imaging of cortical activation. JMRI 5, 49 – 56. Ogawa, S., Tank, D.W., Menon, R., Ellermann, J.M., Kim, S.G., Merkle, H., Ugurbil, K., 1992. Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging. Proc. Natl. Acad. Sci. U. S. A. 89, 5951 – 5955.

D.A. Soltysik et al. / NeuroImage 22 (2004) 1117–1127 Oppenheim, A.V., Schaffer, R.W., Buck, J.R., 1999. Discrete-Time Signal Processing. Prentice Hall, Upper Saddle River, NJ. Peck, K.K., Sunderland, A., Peters, A.M., Butterworth, S., Clark, P., Gowland, P.A., 2001. Cerebral activation during a simple force production task: changes in the time course of the haemodynamic response. NeuroReport 12, 2813 – 2816. Pfeuffer, J., McCullough, J.C., Van de Moortele, P.F., Ugurbil, K., Hu, X., 2003. Spatial dependence of the nonlinear BOLD response at short stimulus duration. NeuroImage 18, 990 – 1000. Rees, G., Howseman, A., Josephs, O., Frith, C.D., Friston, K.J., Frackowiak, R.S.J., Turner, R., 1997. Characterizing the relationship between BOLD contrast and regional cerebral blood flow measurements by varying the stimulus presentation rate. NeuroImage 6, 270 – 278.

1127

Rieke, F., Warland, D., Ruyter van Stevenick, R., Bialek, W., 1997. Spikes: Exploring the Neural Code. MIT Press, Cambridge. MA. Robson, M.D., Dorosz, J.L., Gore, J.C., 1998. Measurements of the temporal fMRI response of the human auditory cortex to trains of tones. NeuroImage 7, 185 – 198. Singh, M., Kim, S., Kim, T.S., 2003. Correlation between BOLD-fMRI and EEG signal changes in response to visual stimulus frequency in humans. Magn. Reson. Med. 49, 108 – 114. Takanashi, M., Abe, K., Yanagihara, T., Oshiro, Y., Watanabe, Y., Tanaka, H., Hirabuki, N., Nakamura, H., Fujita, N., 2000. The effect of stimulus presentation rate on the activity pattern of primary somatosensory cortex; an fMRI study. Proceedings of 8th ISMRM, Denver, p. 996. Vazquez, A.L., Noll, D.C., 1998. Nonlinear aspects of the BOLD response in functional MRI. NeuroImage 7, 108 – 118.