The neural basis of the hemodynamic response nonlinearity in human primary visual cortex: Implications for neurovascular coupling mechanism

Rapid Communication

www.elsevier.com/locate/ynimg NeuroImage 32 (2006) 616 – 625

The neural basis of the hemodynamic response nonlinearity in human primary visual cortex: Implications for neurovascular coupling mechanism Xiaohong Wan,a,b,* Jorge Riera,a Kazuki Iwata,a Makoto Takahashi,b Toshio Wakabayashi,b and Ryuta Kawashimaa a

Advanced Science and Technology of Materials, NICHe, Tohoku University, Sendai, 980-8579 Miyagi, Japan Department of Quantum Science and Energy Engineering, Tohoku University, Sendai, 980-8579 Miyagi, Japan

b

Received 1 September 2005; revised 15 March 2006; accepted 23 March 2006 Available online 11 May 2006 It has been well recognized that the nonlinear hemodynamic responses of the blood oxygenation level-dependent (BOLD) functional MRI (fMRI) are important and ubiquitous in a series of experimental paradigms, especially for the event-related fMRI. Although this phenomenon has been intensively studied and it has been found that the post-capillary venous expansion is an intrinsically nonlinear mechanical process, the existence of an additional neural basis for the nonlinearity has not been clearly shown. In this paper, we assessed the correlation between the electric and vascular indices by performing simultaneous electroencephalography (EEG) and fMRI recordings in humans during a series of visual stimulation (i.e., radial checkerboard). With changes of the visual stimulation frequencies (from 0.5 to 16 Hz) and contrasts (from 1% to 100%), both the event related potentials (ERPs) and hemodynamic responses show nonlinear behaviors. In particular, the mean power of the brain electric sources and the neuronal efficacies (as originally defined in the hemodynamics model [Friston et al. Neuroimage, 12, 466 – 477, 2000], here represent the vascular inputs) in primary visual cortex consistently show a linear correlation for all subjects. This indicates that the hemodynamic response nonlinearity found in this paper primarily reflects the nonlinearity of underlying neural activity. Most importantly, this finding underpins a nonlinear neurovascular coupling. Specifically, it is shown that the transferring function of the neurovascular coupling is likely a power transducer, which integrates the fast dynamics of neural activity into the vascular input of slow hemodynamics. D 2006 Elsevier Inc. All rights reserved. Keywords: EEG; fMRI; Neurovascular coupling; Brain electric source; Hemodynamic response

* Corresponding author. Current address: Labaratory for Cognitive Brain Mapping, Brain Science Institute, RIKEN, Wako, Saitama 351-0198, Japan. Fax: +81 48 462 4651. E-mail address: [email protected] (X. Wan). Available online on ScienceDirect (www.sciencedirect.com). 1053-8119/$ - see front matter D 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2006.03.040

Introduction Since its inception, functional magnetic resonance imaging (fMRI) based on blood oxygenation level-dependent (BOLD) has been extensively used to investigate the sensory, perceptual and recognitive responses in human subjects (Ogawa et al., 1990a,b; Kwong et al., 1992). Generally the BOLD signals are implicitly analyzed based on the assumption that the hemodynamic responses of fMRI are linearly linked with the underlying neural activity (Friston et al., 1994). The evidences which support this assumption have been collected from a variety of experiments in animals and humans (Boynton et al., 1996; Rees et al., 2000; Heeger et al., 2000; Logothetis et al., 2001; Smith et al., 2002). Unfortunately, some data obtained in animal studies have been beyond the linear interpretation (Ogawa et al., 2000; Devor et al., 2003; Jones et al., 2004). These exceptions imply that it is unlikely that a simple linear model is general. If so, we could not always trust BOLD signals to proportionally represent the underlying neural activity. The nonlinear hemodynamic responses using PET (Fox and Raichle, 1984, 1985) and fMRI (Kwong et al., 1992; Boynton et al., 1996; Vazquez and Noll, 1998, Friston et al., 1998; Birn et al., 2001; Singh et al., 2003a,b) in human primary visual cortex with changes of stimulus properties (durations, contrasts and frequencies) have long been studied. Many of the previous studies focused on the BOLD response to different stimulus durations (Boynton et al., 1996; Vazquez and Noll, 1998; Birn et al., 2001). It was found that while longer duration stimuli behave in an approximately linear fashion, short duration stimuli produce responses larger than those predicted from a linear model. These nonlinear effects have been manifested by the fact that the post-capillary venous expansion is an intrinsically nonlinear mechanical process (Buxton et al., 1998; Friston et al., 1998). Nevertheless, the observed fMRI response to a stimulus is the complex result of several cascaded processes, including event evoked neural activity, neurovascular coupling and hemodynamic response (Arthurs and Boniface, 2002). Thus, the hemodynamic response nonlinearity could arise

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

from any other or mixture of these effects, instead of the vascular effect alone. With respect to the hemodynamic response effect, the dynamic correlation between the regional cerebral blood flow (rCBF) and BOLD signals has been extensively investigated by experimental (Zhu et al., 1998; Hoge et al., 1999) and theoretical works (Buxton et al., 1998, 2004; Friston et al., 2000; Friston, 2002; Riera et al., 2004). The Balloon model1 has been suggested to account for the nonlinear convolution of rCBF and BOLD signal (Buxton et al., 1998). Its predictions are consistent with the experimental results measured by the simultaneous BOLD and perfusion imaging technique (Hoge et al., 1999). Since the blood flow is not measurable by BOLD fMRI, Friston et al. (2000) extended the Balloon model by using the stimulus input as an index of neural activity, instead of the cerebral blood flow. However, it should be noted that there are two assumptions underlying this approach: that the neurovascular coupling is linear and that there is a linear relationship between the stimulus input and the event evoked neural activity. As yet, however, the mechanism of neurovascular coupling has been underdetermined. First, it is uncertain whether the hemodynamic response changes in a linear fashion with the underlying neural activity. Earlier studies suggest a predominately linear coupling between BOLD (or CBF) and local field potentials (LFPs) (Logothetis et al., 2001; Lauritzen, 2001), and spiking activity (Rees et al., 2000; Heeger et al., 2000). More recently, the nonlinear coupling of neural activity and hemodynamic response in animals was found by several groups using the optical intrinsic signal imaging (Devor et al., 2003; Jones et al., 2004; Sheth et al., 2004). Nevertheless, in these studies they mixed the neurovascular coupling and hemodynamic response effects together. Second, there is no consensus on which aspect of neural activity drives the hemodynamic response. Logothetis and his colleagues (2001) found the BOLD responses are associated more with the input of an activated region, as reflected by the LFP, than with the spiking output, while others have found a strong correlation also with the spiking activity (Smith et al., 2002; Devor et al., 2005). They argued that the input and output signals are not entirely distinct in a neuronal ensemble. Recently, it has been indicated that the robust correlation of neural activity and the BOLD signal is only at supramillimeter spatial scales (Kim et al., 2004). Although some of the previous studies (Friston et al., 2000; Miller et al., 2001; Birn and Bandettini, 2005) considered both hemodynamic and neuronal nonlinearities, the absence of actual measurements of neuronal activities in these studies makes it difficult to clearly separate the two effects. The advent of the new noninvasive technique with simultaneous electroencephalography (EEG) and fMRI recordings makes it possible to investigate the spatiotemporal correlation between the electrophysiological activities and the hemodynamic responses in human subjects. Although there are many investigations being carried out which take advantage of the concurrent combinational technique, the spatiotemporal correlation between the two modalities has not been well known (Dale and Halgren, 2001; Babiloni et al., 2004, Riera et al., 1 The Balloon model is a hemodynamic input-state-output system model linking between the rCBF (input) and BOLD (output) signals, with two state variables cerebral blood volume (v) and deoxyhemoglobin content ( q). In the context of the Balloon model, the BOLD signal is consisted of an extra and intravascular component, dependent on their deoxyhemoglobin content and weighted by their respective volumes.

617

2005). Partially it is due to that (1) it is difficult to obtain clean event-related potentials (ERPs) from the continuous EEG data contaminated by the artifacts generated by fMRI; (2) the scalp EEG signals suffer from ambiguities in definition of spatial origin of the activities. In this paper, we have used simultaneous high density EEG (64 channels) and 1.5 T fMRI recordings to investigate the spatiotemporal correlation between the BOLD signals and the neural activity in human primary visual cortex with changes of visual stimulation frequencies and contrasts. Our aim is to assess the neural basis of the hemodynamic response nonlinearity and to unveil the mechanism of the neurovascular coupling. An important distinction between the current study and previous studies on the nonlinearity of the BOLD response to different stimulus durations (Boynton et al., 1996; Vazquez and Noll, 1998; Birn et al., 2001) is the type of nonlinearity studied. This study focuses on the BOLD response to different stimulation contrasts and frequencies, with the stimulus duration kept constant. Both types of stimuli (different durations and different contrasts or frequencies) do involve neuronal transients, adaptation, and refractory effects, but the hemodynamic nonlinearities may be different.

Materials and methods Subjects and tasks Eight healthy subjects (three females) initially participated in this study after giving informed consent in accordance with Tohoku University’s Institutional Review Board. Their ages ranged from 22 to 38 years, with a mean age of 27.6 years. Two of them were found to have fallen asleep during experiments and one had unrecoverable EEG signals due to too frequently significant movement artifacts. Therefore, the data from the 5 remaining subjects were examined in the study. There were two types of visual stimulation tasks carried out in this study. A circled radial checkerboard (white and black, 10 rings at 10-eccentricity, checker angle: 10-) was reversed for 4 s (the ‘‘ON’’ period) in each trial (25 s) with changing the stimulation frequencies (0.5, 1.0, 4.0, 8.0 and 16 Hz) under a constant averaged luminance (at the 100% contrast), or with changing the stimulation contrasts (1.0%, 3.2%, 10.0%, 31.6% and 100%, at a log-scale) at the frequency of 2 Hz, followed by a gray screen of equal luminance to the flicking checkerboard for 21 s (the ‘‘OFF’’ period). The subjects were instructed to view the central fixation of a small white filled circle (visual angle: 0.3-) during fMRI scanning, and no overt response was required in any condition. fMRI acquisition The T*2-weighted fMRI images were acquired using a 1.5 Tesla scanner, equipped with a circular head coil (Siemens Magnetom Symphony, Erlangen, Germany), with a gradient echoplanar sequence (repeat time, TR: 2500 ms; echo time, TE: 40 ms; flip angle, FA: 85-; field of view, FOV: 192 mm; voxel size, 3.0  3.0  6.0 mm3; 20 sagittal slices with 0.6 mm gap and bandwidth, 2170 Hz). In both tasks, each subject had three consecutive sessions under five randomly different conditions. Taken together, each subject was presented with a total of 10 trials for each condition. The structural scans were obtained with a spoiled grass gradient sequence with the following parameters: 176 sagittal slices; 1.0 mm thickness with no skip; repeat time,

618

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

TR: 2000 ms; echo time, TE: 4 ms and flip angle, FA: 15-. The images were reconstructed as a 176  256  256 matrix with a 1.0  1.0  1.0 mm3 spatial resolution. EEG recording EEG was simultaneously recorded inside the 1.5 Tesla MRI scanner with a commercially available MR-compatible 64-channel BrainAmp system (Brain Products GmbH, Munich, Germany), sampled at 5 kHz with a band-pass filter (BPF) between 0.5 and 125 Hz. The 56 scalp EEG electrodes were distributed according to the international 10/10 system. The reference was attached to the subjects’ left earlobe. The impedances of all channels were maintained below 5 kV. fMRI data analysis The EEG and fMRI data from each individual subject were analyzed separately. Prior to the correlation analysis, the analysis of the fMRI data was independent of the associated EEG analysis. First, the activated regions were identified and the hemodynamic impulse responses under different conditions were fitted in a conventional analysis (see below). Subsequently, the regional activities obtained from the first analysis were entered into the second Bayesian inference of the hemodynamic system (Friston, 2002). The indices of the BOLD responses derived from the hemodynamic modeling are the vascular inputs, parameterized by the neuronal efficacies (() in this case. The fMRI images preprocessing (realignment, slice timing adjustment, coregistration, spatial normalization and spatial smoothing with 6-mm isotropic full width at half maximum (FWHM) Gaussian kernel) and statistical analysis were carried out using statistical parametric mapping (SPM2, Wellcome Department of Imaging Neuroscience, UCL). Each condition was modeled as a box-car function with a 4-s duration convolved with a canonical, double gamma (one positive and one negative, shifted 2 s apart) hemodynamic response function (HRF) and its time and dispersion derivatives to give three regressors for each condition. The statistical model included global, low frequency and motionrelated confounds and residual serial correlation (Penny et al., 2003). The conjunction analysis of all conditions (Price and Friston, 1997) showed that the activated regions were primarily located in primary visual cortex (Fig. 1, P < 0.001, uncorrected).

Fig. 1. The activated regions (primarily in primary visual cortex) by conjunction analysis of all the conditions (subject 1) with changes of visual stimulation frequencies and contrasts (The color indexes the T value). (For interpretation of the references to colour in this figure legend, the reader is reffered to the web version of this article.)

The dynamics of the activated regions in primary visual cortex was characterized in terms of the first principle component of the adjusted data (removing the global drift and the effects of confounds) from voxels within a 6-mm spherical volume centered on the bilateral maximum of the activated regions. Only voxels surviving a T threshold of P = 0.001 in the first SPM analysis were considered. The stimulation functions modeling the inputs of the hemodynamic modeling comprised a stereotypical train of spikes indexing the presentation of visual stimulation for all conditions (Here we selected the value at the maximum frequency stimulation of 16 Hz reversing at the duration of 4 s, so that the number of spikes is 64). In the hemodynamic modeling (Friston et al., 2000), multiplying of neuronal efficacy and stimulation function is equivalent to the vascular input from the neuronal activity to the ensuing hemodynamics. In Friston’s original application of the hemodynamic modeling to the word listening task at different stimulation rates, the stimulation function u(t) was used as a function of the word presentation rate, and then the estimated neuronal efficacy was assumed to be the same at different presentation rates. However, in our current application, even in the cases of changing stimulation frequencies, the stimulation function u(t) was fixed. Therefore, the neuronal efficacy here is equivalent to the vascular input. The neuronal efficacy was estimated using the Bayesian inference (Friston, 2002, the spm_hdm_ui code in SPM2, we adapted it for multiple sessions). EEG data analysis The imaging and ballistocardiogram artifacts generated by fMRI scanning in the EEG data were removed using the adaptive FIR method (Wan et al., 2006a) and the Wavelet-based Nonlinear Noise Reduction (WNNR) method (Wan et al., 2006b), respectively. Using these two methods, the temporal EEG signals were able to be accurately retrieved. The ERP during the visual stimulation of each condition was calculated by averaging the 10 trials of corrected EEG data with an artifact rejection level of T50Av. The wavelet denoising method (Quian Quiroga and Garcia, 2003) was used to further reduce the ERP noises. A 3-D digitizer (Polhemus Isotrack) was used to detect the electrode positions over the subject’s scalp outside the scanner room. The reconstructed electrodes’ coordinates were matched to the head surface extracted from the subject’s own T1-weighted MRI image. Making a link between the scalp EEG and BOLD signals requires the scalp EEG potentials to be transformed into the brain electric sources in the same space of the BOLD signals. Technically, the fMRI-constrained distributed current source solutions were obtained using the linearly constrained minimum variance (LCMV) method (Van Veen et al., 1997; Robinson and Vrba, 1999; Singh et al., 2003a). The gain matrix of the forward problem was created using the boundary element method (BEM) based on the subject’ realistic head model (Darvas et al., 2004). The source grids (4.0  4.0  4.0 mm3) were generated in the activated regions inferred by SPM2. Subsequently, the mean EEG sources were calculated by averaging the current sources within the regions of primary visual cortex selected in the fMRI analysis. Finally, the mean powers of the estimated EEG current sources during the stimulation period ( T1 ~Tt¼1 jCSDt j2 , T is the stimulation duration, 4000 ms, t is the time point) were used to make a linear regression analysis with the neuronal efficacies inferred from the BOLD data.

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

Results In the present study, the hemodynamic and electrophysiological responses in primary visual cortexes of human subjects were simultaneously recorded by fMRI and EEG with changes of the visual stimulation frequencies and contrasts. Examples of the results for one subject are shown in Fig. 2. The amplitude of the mean event-related hemodynamic response which was calculated by averaging those estimated hemodynamic responses within the selected activated regions (see Materials and methods) increases as the stimulus frequency increases till 8 Hz, and decreases at the frequency of 16 Hz (Fig. 2a). To test whether this changing trend is statistically significant, we compared the five frequency conditions both with each other and with the subgroups. Only the difference between the high (4, 8 and 16 Hz) and low (0.5 and 1.0 Hz) frequency group achieved statistical significance across all five subjects ( P < 0.05, paired t test). In comparison, the mean eventrelated hemodynamic response monotonically increases as the stimulation contrast increases from 1.0% to 100% (Fig. 2b), but increases much more slowly at the high contrast region, presumably due to the saturation effect. Both the nonlinear event-related hemodynamic responses to varying the frequencies and contrasts are in good agreement with previous works (Tootell et al., 1995; Singh et al., 2003b). Simultaneously observed ERPs (the POz channel) are shown in Figs. 2c and d. In Fig. 2c, all EEG responses to changing the frequencies are shown to be synchronized with the reversing visual stimulation (see also Fig. 4). Except for the stimulation frequency of 0.5 Hz, significant attenuation of the second and subsequent EEG responses (N1, 125 ms latency at the transition states) was observed within the 4 s stimulation, most likely reflecting neuronal inhibitory interactions (Simons, 1985). The amplitudes of the late attenuated EEG responses decrease with increasing stimulation frequencies, while the amplitudes of the first EEG responses are very similar across all stimulation frequencies (see integration of the first EEG responses in Fig. 2c). These scalp EEG responses are similar to the LFP observations in rat somatosensory cortex (Sheth et al., 2004). Unlike the conventional analyses of LFP and multiple unit activity (MUA) (Lauritzen, 2001; Sheth et al., 2004), the integrated electrophysiological activity produced by a stimulus was characterized by integration 2 of the ERP power during the stimulation interval (~M m¼1 jERPm j ). In accordance with the hemodynamic responses, a similar integrated electrophysiological response curve against the stimulation frequencies is depicted in the last figure of Fig. 2c. In Fig. 2d, the amplitudes of the ERPs are shown to consistently decrease with decreasing the stimulation contrasts, which also occur with the corresponding hemodynamic responses. This is consistent with another study (Zaletel et al., 2004). It is interesting to note that the ERP latency slightly increases with decreasing the stimulation contrast. Due to the same flickering frequency of 2 Hz used in the contrast conditions, the peak amplitudes of the average ERPs for a spike response are plotted against the stimulation contrasts, as shown in the inset of Fig. 2d. The consistent changes of the hemodynamic and electrophysiological responses to the stimulation frequencies and contrasts suggest that the hemodynamic response nonlinearity probably has its neural origin. However, the hemodynamic responses themselves have nonlinear properties (Buxton and Frank, 1997). In order to separate the two effects, we employed the hemodynamic modeling (Friston et al., 2000; Friston, 2002), which also accounts for the later hemodynamic response nonlinearity, to obtain the neuronal efficacy, representing

619

the vascular input driving the BOLD response. Finally, to quantitatively determine the relationship between the electrophysiological activities and hemodynamic responses, a linear correlation analysis was carried out between the mean powers of the EEG current sources and the neuronal efficacies. The regressive results show that their linear relationship is statistically significant for all 10 conditions (Fig. 3a, R 2 = 0.96, Pearson test) and all 5 subjects (Fig. 3b and Table 1). It is worthy noting that the slopes of the linear regression among the five subjects are clearly different. This may be due to the difference of the hemodynamic response efficiency across the subjects (Logothetis and Wandell, 2004). The extrapolated regression lines consistently cross the y axis with intercepts close to zero.

Discussion Our results show that the hemodynamic response nonlinearity with changes of the visual stimulation frequencies and contrasts is mainly induced by the nonlinear electrophysiological activity, simultaneously measured by EEG, rather than the refractory of hemodynamic responses (Friston et al., 1998). The highly linear correlation between the mean powers of the EEG current sources and the neuronal efficacies derived from fMRI data demonstrates that the fMRI and EEG signals are consistent at the large scale in a variety of conditions. Although the BOLD signals with high spatial resolutions have been linked linearly with the local electrophysiological activity, measured by LFP and/or MUA (Logothetis et al., 2001; Rees et al., 2000), there is recent evidence showing the relationship between the hemodynamic responses and the local electrophysiological activities is nonlinear (Devor et al., 2003; Jones et al., 2004; Sheth et al., 2004). This may be due to the fact that the hemodynamic signal at a given position integrates the electrophysiological activity over a broader spatial region, likely larger than the location with the LFP and MUA (Devor et al., 2005) and probably at supra-millimeter scales (Kim et al., 2004). This renders it promising to correlate the distributed EEG sources with BOLD signals voxel by voxel if the inverse solution of EEG sources with several millimeters spatial resolutions could be perfectly coregistered with the BOLD activated regions. However, EEG generally has a lower spatial resolution than BOLD fMRI and the two maps usually mismatch at the range of several to tens of millimeters (Disbrow et al., 2000; Towle et al., 2003). Hence, it is probably true that the correlation between the distributed EEG sources and BOLD signals might be sensible only at a large scale, where the mean activity was extracted from the regions of interest (ROI), rather than a voxel by voxel. Although the relationship between the scalp EEG and the LFP/MUA is not yet understood well (Nunez, 1995), both two electrophysiological signals recording the extracellular potentials with different scales represent the population postsynaptic activity in the pyramidal cells (Lopes da Silva and Storm van Leeuwan, 1978). The convergence of these signals (EEG, LFP and MUA) with the hemodynamic responses (the current study and Logothetis et al., 2001; Lauritzen, 2001) supports the neural basis of the hemodynamic responses. However, to date, it remains difficult to distinguish which aspect of the electrophysiological activity (i.e., presynaptic signals, likely associated with neurotransmitters release, or postsynaptic signals, associated with neuronal membrane currents, or action potentials; see also below) is better reflected by the hemodynamic responses under normal

620

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

Fig. 2. The hemodynamic responses and electrophysiological activities with changes of stimulation frequencies and contrasts (subject 1): (a) The hemodynamic impulse responses (HIR) (mean T SE) to varying the stimulation frequencies (0.5, 1.0, 4.0, 8.0, 16.0 Hz), averaged within the regions of interest (ROI, 6-mm sphere centered at the maximum of the activated voxels in primary visual cortex). The inset shows the HIR maximum peaks change with the stimulation frequencies; (b) HIR (mean T SE) to varying the stimulation contrasts (1.0%, 3.2%, 10.0%, 31.6%, 100%, at 2 Hz), averaged within the same ROI as in panel a. The inset shows the HIR maximum peaks change against the stimulation contrasts. The simultaneous electrophysiological activity (ERPs, POz) with changes of the stimulation frequencies and contrasts is showed in panels c and d, respectively. The original EEG data were processed to remove the imaging and ballistocardiogram artifacts, and then the ERPs were calculated by averaging across the 10 trials for each condition after artifact rejection. The curves of the last figure in panel c show the integration of the ERP powers with the stimulation frequencies, and the inset in panel d shows the ERP peaks with the stimulation contrasts. The error bars indicate the standard errors.

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

Fig. 2 (continued ).

621

622

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

Fig. 3. The correlation between the mean powers of EEG current sources ( T1 ~Tt¼1 jCSDt j2 ) and the neuronal efficacies (() derived from the BOLD signals in the ROI of primary visual cortex. (a) The neuronal efficacies against the mean powers of EEG current sources across all 10 conditions in subject 1 were well fitted (R 2 = 0.96, P < 0.05, Pearson test) by a linear regression with a slope of 0.143 and an intercept of 0.08. The results of the other subjects are listed in Table 1. (b) The linear correlation between the neuronal efficacies and the mean powers of EEG current sources across all 10 conditions and 5 subjects (R 2 = 0.91, P < 0.05, Pearson test). The error bars indicate the standard errors.

activity (Logothetis and Wandell, 2004). This means the hemodynamic signal not only integrates the neural activity in the spatial region, but also in the temporal interval. This concept is reflected by the observation that the EEG response dynamics sharply varied, but the discrepancy of the BOLD signals slightly changed, especially in the case of varying the stimulation frequencies. Nevertheless, the neurovascular coupling could be a more sophisticated transferring function rather than a simple linear LPF. Otherwise, the broad-banded high-frequency electrophysiological activity input could not be effectively transferred to the narrow-banded low-frequency hemodynamic response output (see Fig. 4 in this paper and Fig. 2 in Logothetis et al., 2001). From our experimental results, the neurovascular coupling might be considered as a function of a power transducer which transfers the power of the neural activity across all frequencies into the vascular input. However, the dynamic feature of the neurovascular coupling remains unknown. For example, we do not know whether it works continuously or discretely (Fig. 5). If it works discretely (a timeintegrated and all-or-none response with a critical threshold, likely associated with neurotransmitter release (Kandel et al., 2000; Devor et al., 2003), then what is the critical time window (DTc) of the electrophysiological integration? Until now, most studies have only focused on the spatial correlation with an arbitrary time window of integration. For example, Logothetis et al. (2001) used

conditions (Lauritzen, 2001; Heeger and Ress, 2002; Devor et al., 2005), and this deserves further study. The BOLD signals used to map brain activation are based on local changes in the blood flow, accompanied by significantly smaller changes in oxygen consumption (Fox and Raichle, 1986; Ogawa et al., 1990a,b). These changes are induced by regional cerebral energy consumption, which is coupled with neural activity (Shulman and Rothman, 1998). Hence, our observation that the neuronal efficacies of the BOLD signals is linearly linked with the mean powers of EEG sources, which are likely proportional to the energy consumption of cellular activity, is in parallel with the above notion. Notably, our results are close to the power-law function fitted between the electrophysiological measurements (LFP/MUA) and the hemodynamic responses (HbO: oxyhemoglobin; HbT: total hemoglobin) by spectroscopic optical imaging (Devor et al., 2003). Actually in the report by Logothetis et al. (2001), a linear correlation analysis was also made between BOLD signals and the mean spectral powers of the LFP/MUA activity averaged across all frequencies. Hence, our finding of the linear relationship between the mean powers of the EEG current sources and the neuronal efficacies of BOLD responses is consistent with their earlier studies. The hemodynamic response of BOLD fMRI is generally believed to be a low-pass filter (LPF) of the underlying neural

Table 1 Activated location of the specified region in primary visual cortex (x, y, and z in millimeters) in the standard space defined by Talairach and Tournoux (1988) and the Bayesian inference of the hemodynamic response in the selected regions Subject

1

2

3

4

5

Region definition Left V1 (T max) Right V1 (T max)

4, 94, 4 (9.93) 8, 94, 4 (8.72)

4, 92, 6 (9.65) 6, 94, 6 (10.32)

12, 92, 2 (7.41) 12, 96, 2 (6.46)

8, 92, 4 (8.47) 10, 94, 4 (9.58)

10, 94, 0 (7.82) 12, 94, 0 (8.45)

Statistical inference Voxels b intercept R 2 ( P < 0.05)

113/111 0.143 T 0.009 0.080 T 0.029 0.96

96/99 0.156 T 0.008 0.065 T 0.021 0.98

54/51 0.122 T 0.011 0.092 T 0.030 0.92

86/89 0.168 T 0.007 0.078 T 0.019 0.97

67/69 0.148 T 0.012 0.089 T 0.027 0.90

V1: primary visual cortex, by virtue of the region definition and spatial smoothing, the cortical region (V1) subsumes portions of V2/V3; the number of voxels (left V1/ Right V1) selected for the Bayesian inference (see Materials and methods), the fitted slope of the linear regression (h) (mean T SE) and the linear correlation coefficient square (R 2) with Pearson test.

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

623

Fig. 4. The time-frequency analysis of the mean EEG current sources averaged in the ROI of primary visual cortex across all subjects, when the visual stimulation frequency is changed from 0.5 to 16 Hz. The EEG responses are synchronized with the flicker reversing. Both in the transition states and steady states of EEG activity, the EEG spectrum has high frequency components, which are usually beyond the region of hemodynamic responses (below 1 Hz). The time-frequency powers were calculated using Morlet wavelet analysis. The wavelets used have a 2-D Gaussian profile, such that the full-width half-maximum (FWHM) of the response in the time and frequency domains are Wt and Wf, respectively. We used a wavelet family parameter of 10, which at each frequency, f, gives Wf = 0.235f and Wt = 3.74/f.

the sampling interval (0.25 s) of fMRI as the time window for the electrophysiological integration, while Sheth et al. (2004) integrated the electrophysiological activity over the time window (2.0 s) of the stimulation interval, similar to the current study. Unfortunately, the fMRI data in this study do not have enough SNR to distinguish the hemodynamic response discrepancy stemming from the difference of the electrophysiological integration. Certainly, using simultaneous high-density EEG and higher field strength fMRI recordings will be helpful to understand the mechanism of neurovascular coupling in more detail with some more elegant

experimental designs. In addition, the bottom-up theoretical modeling is helpful to understand the neurovascular coupling mechanism (Riera et al., in press). Finally, the simple fMRI-constrained distributed inverse solution used in this paper does bias the source localization to the fMRI activated regions. Without fMRI constraints, the free source localization would be more diffusive. However, the absence of fMRI constraints does not imply that the source localization is more accurate due to the non-unique solutions of EEG source inverse problem. How to better use fMRI information to enhance

Fig. 5. The continuous or discrete mechanism of the neurovascular coupling that transfers the power of neural activity into the vascular input.

624

X. Wan et al. / NeuroImage 32 (2006) 616 – 625

the EEG source localization by spatial filtering is currently being developed in our group (Wan and Kawashima, submitted for publication). In summary, we used simultaneous fMRI/EEG recordings to study the correlation of event-evoked hemodynamic responses and electrophysiological activity in human primary visual cortex. Our findings support the convergence of the mean EEG and BOLD signals at a large scale. This study demonstrates the neural basis of the hemodynamic response nonlinearity when changing the visual stimulation frequencies and contrasts and indicates the transferring function of the neurovascular coupling is likely a power transducer, which integrates the fast dynamics of neural activity into the vascular input of slow hemodynamics. However, the dynamic mechanism of the neurovascular coupling remains poorly understood, and we believe further studies of the neurovascular coupling are desirable and will lead us to understand more about the neural basis of fMRI.

Acknowledgments We thank Mr. Yokoyama S. and Ms. Sakai Y. for their technical assistance. This study has been supported by Grant-inAid for Scientific Research (C) No. 15500193, JSPS; JST/ RISTEX, R&D promotion scheme for regional proposals promoted by TAO; and the Tohoku University 21st Century Center of Excellence (COE) Program (Ministry of Education, Culture, Sports, Science and Technology) entitled ‘‘A Strategic Research and Education Center for an Integrated Approach to Language, Brain and Cognition’’. Wan X. acknowledges the support of Tohoku University 21COE Program ‘‘Future Medical Engineering Based on Bio-nanotechnology’’.

References Arthurs, O.J., Boniface, S., 2002. How well do we understand the neural origins of the fMRI BOLD signal? Trends Neurosci. 25 (1), 27 – 31. Babiloni, F., Mattia, D., Babiloni, C., Astolfi, L., Salinari, S., Basilisco, A., Rossini, P.M., Marciani, M.G., Cincotti, F., 2004. Multimodal integration of EEG, MEG and fMRI data for the solution of the neuroimage puzzle. Magn. Reson. Imaging 22 (10), 1471 – 1476. Birn, R.M., Bandettini, P.A., 2005. The effect of stimulus duty cycle and ‘‘off’’ duration on BOLD response linearity. NeuroImage 27, 70 – 82. 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. 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 (13), 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, 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. Buxton, R.B., Uludag, K., Dubowitz, D.J., Liu, T.T., 2004. Modeling the hemodynamic response to brain activation. NeuroImage 23, s220 – s233. Dale, A.M., Halgren, E., 2001. Spatiotemporal mapping of brain activity by integration of multiple imaging modalities. Curr. Opin. Neurobiol. 11 (2), 202 – 208. Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R.M., 2004.

Mapping human brain function with MEG and EEG: methods and validation. NeuroImage 23, S289 – S299. Devor, A., Dunn, A.K., Andermann, M.L., Ulbert, I., Boas, D.A., Dale, A.M., 2003. Coupling of total hemoglobin concentration, oxygenation, and neuronal activity in rat somatosensory cortex. Neuron 39, 353 – 359. Devor, A., Ulbert, I., Dunn, Narayanan, S.N., Jones, S.R., A.K., Andermann, M.L., 2005. Coupling of the cortical hemodynamic response to cortical and thalamic neuronal activity. Proc. Natl. Acda. Sci. U. S. A. 102, 3822 – 3827. Disbrow, E.A., Slutsky, D.A., Roberts, T.P.L., Krubitzer, L.A., 2000. Functional MRI at 1.5 Tesla: a comparison of the blood oxygenation level-dependent signal and electrophysiology. Proc. Natl. Acda. Sci. U. S. A. 97, 9718 – 9723. 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. Fox, P.T., Raichle, M.E., 1985. Stimulus rate determines regional brain blood flow in striate cortex. Ann. Neurol. 17 (3), 303 – 305. Fox, P.T., Raichle, M.E., 1986. Focal physiological uncoupling of cerebral blood flow and oxidative metabolism during somatosensory stimulation in human subjects. Proc. Natl. Acda. Sci. U. S. A. 83, 1140 – 1144. Friston, K.J., 2002. Bayesian estimation of dynamical systems: an application to fMRI. NeuroImage 16, 513 – 530. 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 (1), 41 – 52. Friston, K.J., Mechelli, A., Turner, R., Price, C.J., 2000. Nonlinear responses in fMRI: the Balloon model, Volterra kernels, and other hemodynamics. NeuroImage 12, 466 – 477. Heeger, D.J., Ress, D., 2002. What does fMRI tell us about neuronal activity? Nat. Rev., Neurosci. 3, 142 – 151. Heeger, D.J., Huk, A.C., Geisler, W.S., Albrecht, D.G., 2000. Spikes versus BOLD: what does neuroimaging tell us about neuronal activity? Nat. Neuroscience 3, 631 – 633. Hoge, R.D., Atkinson, J., Gill, B., Crelier, G.R., Marrett, S., Pike, G.B., 1999. Linear coupling between cerebral blood flow and oxygen consumption in activated human cortex. Proc. Natl. Acad. Sci. U. S. A. 96, 9403 – 9408. Jones, M., Hewson-Stoate, N., Martindale, J., Redgrave, P., Mayhew, J., 2004. Nonlinear coupling of neuronal activity and CBF in rodent barrel cortex. NeuroImage 22, 956 – 965. Kandel, E.R., Schwartz, J.H., Jessell, T.M. (Eds.), Principles of Neural Science, 4th edR McGraw-Hill Press, New York, pp. P253 – P279. Kim, D.-S., Ronen, I., Olman, C., Kim, S.-G., Ugurbil, K., Toth, L.J., 2004. Spatial relationship between neuronal activity and BOLD functional MRI. NeuroImage 21, 876 – 885. Kwong, K.K., Belliveau, J.W., Chesler, D.A., Goldberg, I.E., Weiskoff, R.M., Poncelet, B.P., Kennedy, D.N., Hoppel, B.E., Cohen, M.S., Turner, R., et al., 1992. Dynamic magnetic resonance imaging of human brain activity during primary sensory stimulation. Proc. Natl. Acad. Sci. U. S. A. 89, 5675 – 5679. Lauritzen, M., 2001. Relationship of spikes, synaptic activity, and local changes of cerebral blood flow. J. Cereb Blood Flow Metab. 21 (12), 1367 – 1383. Logothetis, N.K., Wandell, B.A., 2004. Interpreting the BOLD signal. Annu. Rev. Physiol. 66, 735 – 769. Logothetis, N.K., Pauls, J., Augath, M., Trinath, T., Oeltermann, A., 2001. Neurophysiological investigation of the basis of the fMRI signal. Nature 412, 150 – 157. Lopes da Silva, F., Storm van Leeuwan, W., 1978. The cortical alpha rhythm in dog: the depth and surface profile of phase. In: Brazier, M.A.B., Petsche, H. (Eds.), Architectonics of the Cerebral Cortex, pp. Raven ress. Miller, K.L., Luh, W.M., Liu, T.T., Martinez, A., Obata, T., Wong, E.C.,

X. Wan et al. / NeuroImage 32 (2006) 616 – 625 Frank, L.R., Buxton, R.B., 2001. Nonlinear temporal dynamics of the cerebral blood flow response. Hum. Brain Mapp. 13 (1), 1 – 12. Nunez, P.L., 1995. Neocortical Dynamics and Human Brain Rhythms. Oxford University Press, Oxford. Ogawa, S., Lee, T.M., Kay, A.R., Tank, D.W., 1990a. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc. Natl. Acad. Sci. U. S. A. 87, 9868 – 9872. Ogawa, S., Lee, T.M., Kay, A.R., Tank, D.W., 1990b. Oxygenationsensitive contrast in magnetic resonance image of rodent brain at high magnetic fields. Magn. Reson. Med. 14, 68 – 78. Ogawa, S., Lee, T.M., Stepnoski, R., Chen, W., Zhu, X., Ugurbil, K., 2000. An approach to probe some neural systems interaction by functional MRI at neural time scale down to milliseconds. Proc. Natl. Acad. Sci. U. S. A. 97, 11026 – 11031. Penny, W.D., Kiebel, S., Friston, K.J., 2003. Variational Bayesian Inference for fMRI time series. NeuroImage 19 (3), 727 – 741. Price, C.J., Friston, K.J., 1997. Cognitive conjunction: a new approach to brain activation experiments. NeuroImage 5, 261 – 270. Quian Quiroga, R., Garcia, H., 2003. Single-trial event-related potentials with wavelet denoising. Clin. Neurophysiol. 114, 376 – 390. Rees, G., Friston, K.J., Koch, C., 2000. A direct quantitative relationship between the functional properties of human and macaque V5. Nat. Neusci. 3, 716 – 723. Riera, J.J., Watanabe, J., Iwata, K., Miura, N., Aubert, E., Ozaki, T., Kawashima, R., 2004. A state-space model of the hemodynamic approach: non-linear filtering of BOLD signal. NeuroImage 21, 547 – 567. Riera, J., Aubert, E., Iwata, K., Kawashima, R., Wan, X., Ozaki, T., 2005. Fusing EEG and fMRI based on a bottom-up model: Inferring activation and effective connectivity in neural masses. Philos. Trans. R. Soc. London, Ser. B 360, 1025 – 1041. Riera, J.J., Wan, X., Jimenez, J.C., Kawashima, R., in press. Nonlinear local electro-vascular coupling part I: a theoretical model. Hum. Brain Mapp. Robinson, S.E., Vrba, J., 1999. Functional neuroimaging by synthetic aperture magnetometry (SAM). Recent Advances in Biomagnetism. Tohoku University Press, Sendai, pp. 302 – 305. Sheth, S.A., Nemoto, M., Guiou, M., Walker, M., Pouratian, N., Toga, A.W., 2004. Linear and nonlinear relationships between neuronal activity, oxygen metabolism, and hemodynamic responses. Neuron 42, 347 – 355. Shulman, R.G., Rothman, D.L., 1998. Interpreting functional imaging studies in terms of neurotransmitter cycling. Proc. Natl. Acda. Sci. U. S. A. 95, 11993 – 11998. Simons, D.J., 1985. Temporal and spatial integration in the rat SI vibrissa cortex. J. Neurophysiol. 54, 615 – 635.

625

Singh, K.D., Barnes, G.R., Hillebrand, A., 2003a. Group imaging of taskrelated changes in cortical synchronisation using nonparametric permutation testing. NeuroImage 19 (4), 1589 – 1601. Singh, M., Kim, S., Kim, T.S., 2003b. Correlation between BOLD-fMRI and EEG signal changes in response to visual stimulus frequency in humans. Magn. Reson. Med. 49, 108 – 114. Smith, A.J., Blumenfeld, H., Behar, K.L., Rothman, D.L., Shulman, R.G., Hyder, F., 2002. Cerebral energetics and spiking frequency: the neurophysiological basis of fMRI. Proc. Natl. Acad. Sci. U. S. A. 99, 10765 – 10770. Talairach, P., Tournoux, J., 1988. A Stereotactic Coplanar Atlas of the Human Brain. Thieme, Stuttgart. Tootell, R.B., Reppas, J.B., Kwong, K.K., Malach, R., Born, R.T., Brady, T.J., Rosen, B.R., Belliveau, J.W., 1995. Functional analysis of human MT and related visual cortical areas using magnetic resonance imaging. J. Neurosci. 15 (4), 3215 – 3230. Towle, V.L., Khorasani, L., Uftring, S., Pelizzari, C., Erickson, R.K., Spire, J.P., Hoffmann, K., Chu, D., Scherg, M., 2003. Noninvasive identification of human central sulcus: a comparison of gyral morphology, functional MRI, dipole localization, and direct cortical mapping. NeuroImage 19, 684 – 697. Van Veen, B.D., Van Drongelen, W., Yuchtman, M., Suzuku, A., 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9), 867 – 880. Vazquez, A.L., Noll, D.C., 1998. Non-linear aspects of the blood oxygenation response in functional MRI. NeuroImage 7, 108 – 118. Wan, X., and Kawashima, R., to be submitted. Optimal resolution of EEG/MEG source imaging by spatial filtering. Wan, X., Iwata, K., Riera, J., Kitamura, M., Kawashima, R., 2006a. Artifact reduction for simultaneous EEG/fMRI recording: adaptive FIR reduction of imaging artifacts. Clin. Neurophysiol. 117 (3), 681 – 692. Wan, X., Iwata, K., Riera, J., Ozaki, T., Kitamura, M., Kawashima, R., 2006b. Artifact reduction for EEG/fMRI recording: nonlinear reduction of ballistocardiogram artifacts. Clin. Neurophysiol. 117 (3), 668 – 680. Zaletel, M., Strucl, M., Rodi, Z., Zvan, B., 2004. The relationship between visually evoked cerebral blood flow velocity responses and visualevoked potentials. NeuroImage 22, 1784 – 1789. Zhu, X.H., Kim, S.G., Andersen, P., Ogawa, S., Ugurbil, K., Chen, W., 1998. Simultaneous oxygenation and perfusion imaging study of functional activity in primary visual cortex at different visual stimulation frequency: quantitative correlation between BOLD and CBF changes. Magn. Reson. Med. 40 (5), 703 – 711.