Interactions between different EEG frequency bands and their effect on alpha–fMRI correlations

NeuroImage 47 (2009) 69–76

Contents lists available at ScienceDirect

NeuroImage j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / y n i m g

Technical Note

Interactions between different EEG frequency bands and their effect on alpha–fMRI correlations J.C. de Munck a,⁎,1, S.I. Gonçalves a,1, R. Mammoliti b, R.M. Heethaar a, F.H. Lopes da Silva c,d a

Brain Imaging Section- Department of Physics and Medical Technology, VU University Medical Centre, De Boelelaan 1117, 1081 HV Amsterdam, The Netherlands Department of Neurological and Cardiovascular Sciences, University of Cagliari, Policlinico Universitario, S.S. 554 km 4.500, 09042 Monserrato (CA), Italy Netherlands Institute for Brain Research, Meibergdreef 33, 1105 AZ, Amsterdam, The Netherlands d Swammerdam Institute for Life Sciences, University of Amsterdam, Amsterdam, The Netherlands b c

a r t i c l e

i n f o

Article history: Received 21 January 2009 Revised 16 March 2009 Accepted 4 April 2009 Available online 17 April 2009

a b s t r a c t In EEG/fMRI correlation studies it is common to consider the fMRI BOLD as filtered version of the EEG alpha power. Here the question is addressed whether other EEG frequency components may affect the correlation between alpha and BOLD. This was done comparing the statistical parametric maps (SPMs) of three different filter models wherein either the free or the standard hemodynamic response functions (HRF) were used in combination with the full spectral bandwidth of the EEG. EEG and fMRI were co-registered in a 30 min resting state condition in 15 healthy young subjects. Power variations in the delta, theta, alpha, beta and gamma bands were extracted from the EEG and used as regressors in a general linear model. Statistical parametric maps (SPMs) were computed using three different filter models, wherein either the free or the standard hemodynamic response functions (HRF) were used in combination with the full spectral bandwidth of the EEG. Results show that the SPMs of different EEG frequency bands, when significant, are very similar to that of the alpha rhythm. This is true in particular for the beta band, despite the fact that the alpha harmonics were discarded. It is shown that inclusion of EEG frequency bands as confounder in the fMRI–alpha correlation model has a large effect on the resulting SPM, in particular when for each frequency band the HRF is extracted from the data. We conclude that power fluctuations of different EEG frequency bands are mutually highly correlated, and that a multi frequency model is required to extract the SPM of the frequency of interest from EEG/fMRI data. When no constraints are put on the shapes of the HRFs of the nuisance frequencies, the correlation model looses so much statistical power that no correlations can be detected. © 2009 Elsevier Inc. All rights reserved.

Introduction In many studies where neurocognitive processes are investigated using fMRI and BOLD signals (e.g. Greicius et al., 2003; Damoiseaux et al., 2006), the so called awake “resting state” is used as the state of reference, although in many of these studies this “resting state” is not objectively controlled. This shortcoming has been compensated using simultaneously recorded EEG signals, assuming that in the awake and resting subject the EEG displays prominent alpha activity. This issue was investigated previously (Goldman et al., 2002, Laufs et al., 2003b, 2006; Moosmann et al., 2003, Gonçalves et al., 2006). Since these studies were aimed at the localisation of the generators of the alpha band, temporal variations in alpha power were convolved with the standard hemodynamic response function, and correlated with the ⁎ Corresponding author. Fax: +31 204444147. E-mail address: [email protected] (J.C. de Munck). 1 First and second authors contributed equally to this manuscript. 1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2009.04.029

fMRI BOLD signals to obtain statistical parametric maps (SPMs) in which the significant alpha generators were indicated. Results of the above mentioned studies showed that generally the alpha rhythm was negatively correlated to the BOLD signal in cortical regions, including e.g. the visual cortex or the pre- and post-central gyrus, whereas positive correlations were found in the thalamus. Large inter-subject variability was found in alpha-BOLD SPMs (Gonçalves et al., 2006, De Munck et al., 2007) which appeared to be mainly related to individual variations in the EEG signal (Gonçalves et al., 2008). The few studies that were performed on other frequency bands (Laufs et al., 2003a; Tyvaert et al., 2008) were based on group analysis and showed contradicting results. In De Munck et al. (2007, 2008), the EEG/fMRI correlation model was extended by extracting the hemodynamic response of the alpha from the data using the GLM, instead of assuming the standard HRF. The impulse response function so obtained was called alpha response function (ARF) to stress that it reflects the response to alpha rhythm observed with the surface EEG instead of a true hemodynamic

70

J.C. de Munck et al. / NeuroImage 47 (2009) 69–76

response to the local electrical input. These ARFs were found to vary over the cortex (De Munck et al., 2007) and the peak time delay was shorter for the thalamus with respect to the cortex. These differences of ARFs might be due to differences in regional vasculature, which is a characteristic that has been noted by Harrison et al. (2002), but may also be due to differences in local electric activities caused by the ignorance of other frequency bands than the alpha band in the correlation model. In the present paper, we consider the EEG/fMRI correlation models used so far as incomplete because they are based on the selection of a single EEG frequency band from which a regressor of interest is derived. Ignoring the other frequency bands, reflecting related or unrelated associated cognitive processes, and their interactions with the frequency of interest, may lead to spurious correlations or to the observed spatial variations of the ARF. Tyvaert et al. (2008) have demonstrated, using the EEG/fMRI resting state data of a few subjects with epilepsy and using the standard HRF model, that accounting for the other frequencies may indeed have an effect on the resulting SPM. The main goal of the present study is to investigate the interaction of different EEG frequency bands and their effects on the SPM systematically by comparing different correlation models with various numbers of degrees of freedom. Contrary to Tyvaert et al. (2008), EEG/fMRI data of healthy subjects are used that were all measured with the same protocol. Basically two models were used: (a) one that we call the Single Frequency Model, where a single frequency of interest is selected and the other frequencies are ignored and (b) another one that we call the Multi Frequency Model and where one frequency band is taken as regressor of interest and all the others are added to the confounders. Furthermore in model (a) we used a Free Hemodynamic Response Function (HRF); in model (b) we used two variations of the model. In model (b1) the Free HRF was used for each frequency component and in model (b2) we used the standard or canonical HRF. Related secondary goals of the present study were to explore the interactions between EEG regressors derived from different frequency bands and electrode sides and to study the hemodynamic response functions and SPMs for other frequency bands than the alpha band. Materials and methods Subjects Co-registered EEG–fMRI data were acquired from 16 healthy subjects (7 males, mean age 27, ± 9 years) while they lay rested in the scanner avoiding to fall asleep. These data were the same as presented in (De Munck et al., 2007). Furthermore, an additional data set was recorded for subject 6.

Acquisition of EEG data The EEG was acquired using an MR compatible EEG amplifier (SD MRI 64, MicroMed, Treviso, Italy) and a cap providing 64 Ag/AgCl electrodes positioned according to the extended 10–20 system. The reference electrode was positioned between Pz and POz. Later, while processing the data, all channels were re-referenced to average reference. Data were acquired at a rate of 2048 Hz and an anti-aliasing hardware filter at 537.6 Hz was used. For the subject's comfort, a pillow of viscoelastic foam was placed under the head. Acquisition of fMRI data For each subject, functional and anatomical images were acquired on a 1.5 T MR scanner (Magnetom Sonata, Siemens, Erlangen, Germany) according to the protocol described in (De Munck et al., 2007). In summary, for the fMRI a T2⁎ weighted EPI sequence was used (TR = 3000 ms, TE = 60 ms, 64 × 64 matrix, FOV = 211 × 211 mm, slice thickness = 3 mm (10% gap), voxel size = 3.3 × 3.3 × 3 mm3) with 24 transversal slices covering the complete occipital, parietal and frontal lobes. During the experiment, 600 volumes were recorded, resulting in 30 min of fMRI data for each subject. The anatomical scan was acquired using a T1-weigthed sequence (MPRAGE, TR = 2700 ms, T1 = 950 ms, TE = 5.18 ms, 256 × 192 × 160 matrix, FOV = 256 × 192 × 240 mm, voxel size = 1.0 × 1.0 × 1.5 mm3) consisting of 160 coronal slices of 1.5 mm thickness. Analysis of EEG data The EEG data were corrected for gradient and pulse artefacts using an algorithm developed in-house and described in Gonçalves et al. (2007). Briefly, the dead time before the acquisition of the first slice as well as the time corresponding to the acquisition of one slice was estimated from the MR gradient artefacts in the EEG. These parameters were then used to build a slice artefact template which was subtracted from the data. This was followed by the subtraction of a volume artefact template. Before subtraction, data were shifted with a sub-sample accuracy using sinc interpolation. The pulse artefact was removed by first applying a clustering algorithm to account for its beat-to-beat variability and by afterwards subtracting from each of the detected pulse artefact events the average artefact template corresponding to the appropriate cluster. Data were subdivided into 600 trials of 3 s each, corresponding to each of the recorded fMRI volumes. Data were re-referenced to average reference. Spectrograms were computed for each channel using the FFT algorithm described in Frigo and Johnson (2005). An additional spectrogram corresponding to the occipital leads was derived by averaging the spectrograms corresponding to O1, O2, POz,

Table 1 Summary of the models used in EEG/fMRI correlation analysis. Model

Columns of R

Columns of S

NEP

Single Frequency Model FreeHRF

10 time shifted versions of the power time variations in any of the EEG bands.

35

Multi Frequency Model FreeHRF

8 time shifted versions of the power time variations in any of the EEG bands.

Multi Frequency Model StdHRF

1 power time variations in any of the EEG bands convolved with the standard HRF.

• • • • • • • • • • • • • •

10 time shifted versions of the HBI regressor 6 RETROICOR regressors 6 motion regressors 3 polynomial trends 8 × 4 time shifted versions of power time variations of all four EEG bands not used as regressor 10 time shifted versions of the HBI regressor 6 RETROICOR regressors 6 motion regressors 3 polynomial trends 4 power time variations of all EEG bands not used as regressor convolved with the standard HRF 10 time shifted versions of the HBI regressor 6 RETROICOR regressors 6 motion regressors 3 polynomial trends

65

30

The entries Columns of R and Columns of S refer to the functions which are used as columns of matrix R and S in Eq. (2) respectively. The entry NEP refers to the total number of estimated parameters.

J.C. de Munck et al. / NeuroImage 47 (2009) 69–76

71

Fig. 1. SPMs of subject 8 for several frequency bands, obtained with the Single Frequency Model (a) and with FDR = 0.01. At this significance level, correlation values range from 0.21 to 0.40. Increasing correlation values are colour-coded from dark red to white. The theta band hardly showed correlation in this subject.

PO3, PO4, PO7, PO8, and two additional electrodes positioned between P5 and PO3, and between P6 and PO4, respectively. The EEG bands were defined as delta (0.1–4 Hz), theta (4.5–8 Hz), alpha (8.5–12 Hz), beta (12.5–36 Hz) and gamma (36.5–100 Hz). Furthermore, the alpha rhythm frequencies were determined individually for each subject based on the maximum (within the alpha band) of the average spectral profile and the alpha rhythm harmonics were removed from the beta and gamma bands. Outliers in spectral time series xn of each band were skipped in a semi-automatic way, based on the determination of the median of xn. Analysis of fMRI data In order to register the EPI to the anatomical scan, the first EPI scan was automatically matched to the anatomical scan using the algorithm of Maes et al. (1997). EPI images were afterwards motion corrected and the motion parameters (in a total of six motion time series) were later used as confounders in the correlation analysis (Friston et al., 1996) between the BOLD signal and the power variations of the alpha rhythm. Motion corrected images were then smoothed using a Gaussian kernel with a standard deviation parameter of 5 mm. Additional confounders (six time series in total) were derived from the pulse artefact markers using the RETROICOR method described in Glover et al. (2000). These confounders accounted for the effect of the heart beat phase in the BOLD signal due to the pulsation of veins, arteries and cavities filled with cerebrospinal fluid. In addition to the phase of the heart beat, the heart rate was also added as confounder, similar to De Munck et al. (2008). This extra confounder was derived from the heart beat intervals (HBI), defined as the time interval between two heart beats as given by the pulse or ECG signals. Due to the fact that a fast response is expected to variations in the HBI, a separate HBI confounder was taken for each slice where the acquisition time of each slice was taken into account. Details over the precise construction of this confounder can be found in De Munck et al. (2008). Finally, constant, linear and quadratic trends were also used as confounders. Formally, the BOLD signal is described according to d = Rθ + Sβ + e

regression coefficients or nuisance parameters. Finally, ɛ is an N × 1 column vector modelling Gaussian noise. In the sequel, the statistical significance of OLS estimated θ̂ will be determined using a partial F-test. The correction for multiple comparisons is obtained by controlling the false detection rate (FDR) (Benjamini and Hochberg, 1995; Genovese et al., 2002). In this paper we investigate three different models derived from Eq. (1) in order to explore the relation between BOLD and EEG derived regressors. These models are summarized in Table 1. (a) In the first model (Single Frequency Model) the formalism is similar to that described in De Munck et al. (2007, 2008). In this case, the regressors are shifted versions of the power time variations in any of the EEG bands, taken separately and the time shifts span from − 2 to 7 samples. When one frequency band is studied, the others are ignored. The confounders are the shifted versions of the HBI, the motion parameters as well as the quadratic, linear and constant trends. In the Single Frequency Model, a total of 35 regression parameters are estimated (K = 10 regressors of interest and L = 25 confounders). (b1) The second model (Multi Frequency Model, FreeHRF) considered is similar to model (a) but in this case, one EEG band is taken as regressor of interest and all the other frequency bands are added to the confounders. To reduce the number of regressors the time shifts span from 0 to 7 samples, i.e. anti-causal activations are ignored and therefore a total of 65 regressor parameters are estimated. (This reduction of regressors is justified by a case study presented in

ð1Þ

where d is the N × 1 column vector containing the BOLD time series with N the number of time points, R is the N × K matrix containing the K regressors of interest as columns, θ is the K × 1 column vector containing the response function or the weight of the frequency band of interest. Furthermore, S is N × L matrix containing the L confounders in columns, and β is the L × 1 column vector containing the associated

Fig. 2. The Fraction of Overlapping Points (FOP) between the SPM of a frequency band (delta, theta, beta and gamma) and the SPM of the alpha band is plotted the Fraction of Activated Points (FAP) of the band under consideration. In this plot an FDR threshold of 0.01 was used and all SPMs were computed using the Single Frequency Model, with FreeHRF, model (a).

72 J.C. de Munck et al. / NeuroImage 47 (2009) 69–76 Fig. 3. The estimated HRFs at the Occ (left column), Lpar (middle column) and Thal (right column) points, detected at an FDR of at most 0.05, are plotted for all frequency bands and shown in different rows. The colours with which the HRFs are plotted correspond to different subjects.

J.C. de Munck et al. / NeuroImage 47 (2009) 69–76

73

order to compare SPMs obtained using different models or frequency bands, several quantities were defined. The fraction of activated points (FAP) of a given SPM is defined as FDR

FAPk =

Nk ; Total NBrain

ð3Þ

is the number of activated points within the brain for where NFDR k model k at a given FDR and NTotal Brain is the total number of brain points in the fMRI scan. In addition, in order to evaluate the similarity between two SPMs obtained with different models k and l, the fraction of overlapping points (FOP) was determined according to: FOPk;l =

FDR Nk;l Total NBrain

;

ð4Þ

is the number of points activated in the SPMs of both where NFDR k models k and l for a given FDR. Note that by definition, one has FOPl,k = FOPk,l ≤ FAPk. In the sequel the model subscripts k and l will be omitted. Analysis of EEG derived regressors Fig. 4. Plot of the estimated HRFs (FDR b 0.05) determined from EEG/fMRI data of subject 6, recorded with time interval of about three months (1st meas and 2nd meas). (A) delta band. (B) theta band.

De Munck et al. (2009) where anti-causal alpha-BOLD correlations disappeared when other frequency bands were added as confounder). Finally, (b2) the third model considered (Multi Frequency Model, StdHRF) uses the standard HRF instead of the free HRF. In this model, all bands not considered as regressor of interest were added as confounders. In the Multi Frequency, StdHRF model 30 regression parameters are estimated and the following HRF model was adopted: 8 tV0 <0

b b1 2 ; ð2Þ hðt Þ = t t : c1 expð−t = t1 Þ + c2 expð−t = t2 Þ t N 0 t1 t2 with t1 = 5 s, t2 = 10 s, b1 = 0.8, b2 = 0.9, c1 = 1 and c2 = −0.3. The models (Table 1) are used to determine the correlation between one of the EEG rhythms and the BOLD time series at each voxel thus generating SPMs containing the significant correlations. In

Spectral profiles corresponding to different frequency bands were used as regressors in the GLM in order to compute the SPMs representing the correlation between the spontaneous EEG power variations and the variations in the BOLD signal. The characteristics of the EEG derived regressors and their relation to the resulting SPMs using different models were studied. The spatial variation of the spectral profiles corresponding to each band was explored. For that, the cross-correlations between the spectral profiles measured in each channel i and the spectral profile determined by the averaging the occipital profiles were computed according to:   xTi xocc corri = qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  T  T ffi xi xi xocc xocc

ð5Þ

where xi is the spectral profile at channel i (with offset removed), and xocc the spectral profile at the same frequency band, but averaged over all occipital channels (and also offset removed). The normalized histogram of the correlation values corresponding to all subjects was computed with a bin size of 0.1.

Fig. 5. Normalized distribution of the cross-correlation values computed between the spectral profiles of the signal recorded at each channel and the average profile derived from the occipital channels, according to Eq. (4), pooled over all subjects.

74

J.C. de Munck et al. / NeuroImage 47 (2009) 69–76

Table 2 The correlations between different occipital frequency band regressors are shown: De = delta, Th = theta, Al = alpha, Be = beta and Ga = gamma. Subject

De–Th

De–Al

De–Be

De–Ga

Th–Al

Th–Be

Th–Ga

Al–Be

Al–Ga

Be–Ga

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Average

0.70 0.42 0.31 0.39 0.20 0.65 0.30 0.64 0.36 0.50 0.52 0.46 0.41 0.40 0.67 0.67 0.48

0.64 0.03 0.22 0.15 −0.36 0.48 0.04 0.26 0.05 0.18 0.03 0.15 0.29 0.17 0.27 0.26 0.18

0.55 0.14 0.21 0.12 −0.11 0.48 0.24 0.17 0.23 0.16 0.23 0.22 0.16 0.23 0.27 0.29 0.22

0.64 0.05 0.06 −0.04 0.31 0.35 0.08 0.55 − 0.01 0.49 0.41 0.30 0.42 0.37 0.31 0.42 0.30

0.78 0.29 0.74 0.41 0.16 0.59 0.50 0.33 0.21 0.23 0.36 0.52 0.49 0.37 0.57 0.27 0.43

0.71 0.69 0.66 0.51 0.28 0.44 0.73 0.52 0.48 0.51 0.52 0.49 0.43 0.59 0.50 0.48 0.53

0.67 0.40 0.15 0.18 0.12 0.28 0.23 0.55 −0.03 0.45 0.27 0.28 0.29 0.40 0.51 0.46 0.33

0.79 0.39 0.69 0.61 0.52 0.51 0.68 0.39 0.43 0.49 0.52 0.53 0.43 0.46 0.62 0.25 0.52

0.66 0.06 0.25 0.31 − 0.11 0.18 0.23 0.34 − 0.02 0.34 0.18 0.04 0.31 0.33 0.36 0.19 0.23

0.63 0.56 0.43 0.30 0.03 0.19 0.29 0.23 0.01 0.32 0.32 0.22 0.18 0.45 0.24 0.23 0.29

In particular, there is a strong correlation between the alpha and the beta band, despite the fact that the alpha harmonics are removed from the beta band.

Selection of points of interest In order to explore the spatial variation of the HRF, four anatomical points of interest in the brain were selected for each subject. These points were located in the visual cortex in the mid-sagittal plane (Occ point), at the left and right central sulcus (LPar and RPar points respectively), and at the thalamus (Thal point). Estimated HRFs at each of these points were averaged over a spherical region with a radius of 5 mm. More detailed information about the location of these points can be found in De Munck et al. (2007). Results First we addressed the question whether the different EEG frequency bands are correlated significantly, or not, with the simultaneously recorded BOLD signal. This was done using the Single Frequency Model (a), where each EEG band was separately correlated to the BOLD signal without using the other bands as confounders. The SPMs were generated for the five different bands and a comparison was made between the SPM of the alpha band and the SPMs of each of the other bands. Typical SPMs are shown in Fig. 1 for subject 8 for whom only the theta band (not shown) did not yield a significant SPM at an FDR = 0.01. It can be seen that when the SPM corresponding to a given frequency was significant, it greatly overlapped with the alpha SPM. The corresponding information is presented for all subjects in a summarized form in Fig. 2 where overlaps of the SPMs of all subjects

and all frequency bands are quantitatively compared with the SPM of the alpha band. In this figure the fraction of overlapping points (FOP) is plotted against the fraction of activated points (FAP). Since the alpha SPM is being used as reference, results are represented only for the remaining bands. These results were obtained at FDR = 0.01 and they show that for the SPMs of all frequency bands, when the FAP of that band increases, also the FOP with the alpha SPM increases. This relation is strongly linear, except for very small FAP values. This means that for all frequency bands other than alpha, when there are activated points in the corresponding SPMs, their locations substantially overlap with those of the alpha band SPM. Second, we examined whether different EEG frequency bands showed similar HRFs in their correlation with BOLD. This was done by plotting those functions (normalized) corresponding to the Occ, LPar and Thal points described in the “Materials and methods” section. Only HRFs that were significant at an FDR = 0.05 were plotted. The HRFs corresponding to the RPar points are not shown since they were very similar to those of the LPar. Results are presented in Fig 3 and show that the shapes of the HRFs of the alpha and beta bands over the cortical points (Occ and LPar) were rather similar and vary relatively little over subjects. The HRFs of the delta, theta and gamma bands differed from those of the alpha and beta bands in that they showed much more variability over subjects, even though the number of subjects for whom these functions were significant at the selected threshold was much less than that for the alpha and beta bands. The variability over subjects occurs in terms of the peak response time as well as in terms of the polarity of the maximum response. The only

Fig. 6. Plot of the FAP for all subjects and for three different models: Single Frequency Model (a) and Multi Frequency Model with Free HRF (b1) or standard HRF (b2).

J.C. de Munck et al. / NeuroImage 47 (2009) 69–76

75

Fig. 7. Alpha SPMs obtained for subject 7 using the Multi Frequency Model with either FreeHRF (A) or StdHRF (B). In the left figure correlation values range from 0.25 to 0.37 and increasing correlation values are colour-coded from dark red to white. The presented correlations correspond to an FDR of 0.01. In the right figure significant correlation values range from −0.20 to −0.15 and from 0.15 to 0.20, corresponding to an FDR of 0.01. The colour scale is presented at the right of panel B. Note that when the free HRF is used, multiple parameters are tested and the SPM shows the combined positive correlation coefficient for all estimated parameters of interest. When a fixed HRF is used, there is a single parameter of interest, of which the sign is absorbed in the correlation coefficient.

subject that shows parietal activity in the theta band is subject 6 and this is also the subject where the occipital correlation in the beta and gamma bands is positive instead of negative. The HRFs of the Thalamus were quite different when compared to those of the cortical points, for all frequency bands. Furthermore, the thalamic HRFs of the alpha and beta bands clearly differed from the corresponding HRFs of the cortical points in the sense that the thalamic HRFs were mostly positive or biphasic. They appeared to peak earlier, at approximately 5 s thus confirming what was also found by De Munck et al. (2007) for the alpha band. Also in the case of the delta, theta and gamma bands, the thalamic HRFs appeared to be mostly positive or biphasic. For subject 6, two data sets were recorded with about three months in between in order to verify the reproducibility of these results. Fig. 4 shows that the reproducibility of the HRFs of the Occ point for the delta and theta bands is rather high. In order to exclude the possibility that the variations in the HRFs are biased by the fact that the regressors are derived from the subset of occipital leads, the spatial variation of the regressors was investigated. This was done by comparing the histograms of the cross-correlations, computed using Eq. (5), between the spectral profiles of the signal recorded at each channel and that derived from the occipital channels. Results are presented in Fig. 5 where it can be seen that for the theta, alpha and beta bands, the cross-correlation distribution peaks at the larger values. This implies that within these bands and for all subjects, but much less for the delta band, the regressors derived from the spectral profiles measured at different locations on the surface of the head and those derived from spectral profiles measured at the occipital region are largely correlated and therefore the subset of electrodes to derive the EEG regressors can be restricted to the occipital electrodes without too much loss of generalization. The correlations between regressors derived from occipital leads in different frequency bands are presented in Table 2. Substantial correlations appear to be present between all frequency bands and between the combinations beta–theta and beta–alpha in particular. Note that the latter correlation cannot be caused by the second and third alpha harmonics that are in the usual definitions part of the beta

band. In our definition the harmonics of the alpha band have been omitted in the beta band. Third, we examined whether the choice of the correlation model had influence on the resulting SPM. The alpha SPMs were computed using both versions of the Multi Frequency Models FreeHRF (b1) and StdHRF (b2). The performance of these models was also compared with that of the Single Frequency Model (a). The FAP of the corresponding alpha SPMs (thresholded at an FDR = 0.01) is presented for the various subjects in Fig. 6. It can be seen that the FAP, obtained using multi frequency model (b1) strongly decreases when compared to that obtained using the model (a); this happens in 13 of the 16 subjects. This shows that when all bands not used as regressor of interest, are added as confounders there is a strong decrease in the FAP for almost all subjects. However, if instead of using the free HRF, the standard HRF is used while all bands not used as regressor of interest are added as confounders, model (b2) the FAP increases in 11 subjects while it hardly changes in the remaining subjects. This is

Fig. 8. The two multi frequency models are compared by plotting the FOP of the SPM generated by the FreeHRF against the FAP of the SPM obtained with the StdHRF. The FDR was set at 0.01.

76

J.C. de Munck et al. / NeuroImage 47 (2009) 69–76

illustrated in Fig. 7 where the alpha SPMs of subject 7 obtained with the two versions of the Multi Frequency Model, (FreeHRF and StdHRF) are shown. It appears that when the Free HRF model is used, the EEG/ fMRI correlation at the parietal lobe is lost, compared to the standard HRF. In order to investigate the similarity between the SPMs generated by models (b2) and (b1) more generally, the FAP obtained with the model (b2) is plotted against the FOP obtained (b1), using an FDR = 0.01. Results are presented in Fig. 8. This figure shows first of all that there is a large inter-subject variability between the FAP of different subjects, when model (b2) is used. Generally, when the FAP is large, the overlap of the SPM obtained with model (b1) also tends to be large, indicating that although the SPM of model (b1) contains less activated voxels (Fig. 6), the map itself is generally a subset of the one obtained with model (b2). Discussion and conclusions Since the influential paper of Goldman et al. (2002), localisation of the alpha band generators using EEG/fMRI has been considered as a breakthrough brain imaging because contrary to MEG/EEG inverse modelling, no assumptions on the number of dipoles, their geometry and spatio-temporal smoothness were needed. The mapping of spatially complex patterns is achieved with a systematic k-space sampling of the MRI scanner, whereas for MEG/EEG spatial blurring from source to sensor space poses a limit on the number of spatial parameters that can be sampled. However, considering that the impulse response function to alpha power variations might vary over the brain (De Munck et al., 2007) and that heart beat variations (Shmueli et al., 2007) and respiration (Birn et al., 2006) are correlated to both alpha power and the fMRI signal, more complicated correlation models are required, with already many more (nuisance) parameters in the GLM. Since the time duration of the EEG/fMRI experiments was still sufficiently long, the statistical power of the correlation model remained sufficient to detect the main characteristics of resting state alpha rhythm (De Munck et al., 2008). In the present paper, interaction of alpha band and the other EEG bands is considered and now the limits of statistical power come at sight. It is shown that, when all frequency bands are modelled with an a priori unknown (but causal) HRF, the SPM of the alpha rhythm may reduce almost completely (see e.g. subjects 5 and 8, Fig. 6) or substantial parts may disappear (Fig. 7). Therefore, despite the fact that the shape of HRF may contain physiological relevant information, one has to reduce the number of free parameters of the model, for instance by adopting the standard HRF for all frequency bands, as we did in model (b2). Remarkable findings of our study are that the SPM and the hemodynamic response function of, in particular, the alpha and beta bands are very similar, despite that fact that in our beta band regressors the harmonics of the alpha band were omitted. These findings follow from Figs. 2 and 3 and Table 2. The implication of this finding is that the oscillators of the alpha and beta bands are switched on and off “simultaneously”, i.e. simultaneously on a time scale of several seconds. Fig. 4 furthermore shows that this finding does not depend on the choice of the occipital electrodes because in particular the beta band variations are highly correlated over all electrodes. Our conclusion that the SPM of the beta band greatly overlaps with that of the alpha band contradicts with the finding of Laufs et al. (2003a). These authors found, using a group analysis, that both beta and theta band correlates to areas such as those belonging to the Default Mode Network (Raichle et al., 2001; Greicius et al., 2003; Raichle and Snyder, 2007). For the theta band, we only found EEG/BOLD correlation in two subjects (Fig. 2), and also these SPMs greatly overlap with the SPM of the alpha band. Similar to Tyvaert et al. (2008), who also adopted a multi frequency model, we find large inter-subject variability in the SPMs of EEG/fMRI

correlation analysis. However, our data shows (e.g. Fig. 4) that within subject reproducibility is high. Therefore, we suggest that the multi frequency model using the canonical HRF is still not a complete model and that more advanced models, giving a better mathematical description of the underlying physiology, are needed to bring results from different subjects in mutual agreement. An important step forward would be to realize that the true unknowns in the EEG/fMRI data sets are not so much the θ-variables in Eq. (1), but rather the electric neural activity in different frequency bands, which reveals itself both in the observed EEG and fMRI signals. These more advanced models would truly integrate EEG and fMRI but they would also bring back the EEG inverse modelling issues into the picture. References Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statist. Soc. B 57 (1), 289–300. Birn, R.M., Diamond, J.B., Smith, M.A., Bandettini, P.A., 2006. Separating respiratoryvariation-related fluctuations from neuronal-activity-related fluctuations in fMRI. NeuroImage 31, 1536–1548. Damoiseaux, J.S., Rombouts, S.A.R.B., Barkhof, F., Scheltens, P., Stam CJ Smith, S.M., Beackmann, C.F., 2006. Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. 103 (37), 13848–13853. De Munck, J.C., Gonçalves, S.I., Huijboom, L., Kuijer, J.P.A., Pouwels, P.J.W., Heethaar, R.M., Lopes da Silva, F.H., 2007. The hemodynamic response of the alpha rhythm: an EEG/ fMRI study. NeuroImage 35, 1142–1151. De Munck, J.C., Gonçalves, S.I., Faes, T.h.J.C., Kuijer, J.P.A., Pouwels, P.J.W., Heethaar, R.M., Lopes da Silva, F.H., 2008. A study of the brain's resting state based on alpha band power, heart rate and fMRI. NeuroImage 42, 112–121. De Munck, J.C., Goncalves, S.I., Van Houdt, P., Ossenblok, P., Lopes da Silva, F.H., 2009. Estimating the hemodynamic response of EEG-features. In: Ullsperger, M., Debener, S. (Eds.), Simultaneous EEG and fMRI: Recording, Analysis and Application. Oxford University Press. Frigo, M., Johnson, S.G., 2005. The design and implementation of FFTW. Proc. of the IEEE 93 (2), 216–231. Friston, K.J., Williams, S., Howard, R., Frackowiak, R.S., Turner, R., 1996. Movementrelated effects in fMRI time-series. Magn. Reson. Med. 35, 346–355. Genovese, C.R., Lazar, N.A., Nichols, T., 2002. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. NeuroImage 15, 870–878. Glover, G.H., Li, T.Q., Ress, D., 2000. Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn. Reson. Med. 44, 162–167. Goldman, R.I., Stern, J.M., Engel, J., Cohen, M., 2002. Simultaneous EEG and fMRI of the alpha rhythm. NeuroReport 13 (18), 2487–2492. Gonçalves, S.I., de Munck, J.C., Pouwels, P.J.W., Schoonhoven, R., Kuijer, J.P.A., Maurits, N.M., Hoogduin, J.M., Van Someren, E.J.W., Heethaar, R.M., Lopes da Silva, F.H., 2006. Correlating the alpha rhythm to BOLD using simultaneous EEG/fMRI: inter-subject variability. NeuroImage 30 (1), 203–213. Gonçalves, S.I., Pouwels, P.J.W., Kuijer, J.P.A., Heethaar, R.M., De Munck, J.C., 2007. Artefact removal in co-registered EEG–fMRI by selective average subtraction. Clin. Neurophysiol. 118 (11), 2437–2450. Gonçalves, S.I., Bijma, F., Pouwels, P.W.J., Jonker, M.A., Kuijer, J.P.A., Heethaar, R.M., Lopes da Silva, F.H., De Munck, J.C., 2008. A data and model-driven approach to explore inter-subject variability of resting-state brain activity using EEG-fMRI. IEEE J. Sel. Top. SP 2 (6), 944–953. Greicius, M.D., Krasnow, B., Reiss, A.L., Menon, V., 2003. Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc. Natl. Acad. Sci. 100 (1), 253–258. Harrison, R.V., Harel, N., Panesar, J., Mount, R.J., 2002. Blood capillary distribution correlates with hemodynamic-based functional imaging in cerebral cortex. Cereb. Cortex 12, 225–233. Laufs, H., Krakow, K., Sterzer, P., Eger, E., Beyerle, A., Salek-Haddadi, A., Kleinschmidt, A., 2003a. Electroencephalographic signatures of attentional and cognitive default modes in spontaneous brain fluctuations at rest. Proc. Natl. Acad. Sci. 100 (19), 11053–11058. Laufs, H., Kleinschmidt, A., Beyerle, A., Eger, E., Salek-Haddadi, A., Preibisch, C., Krakow, K., 2003b. EEG-correlated fMRI of human alpha activity. NeuroImage 19, 1463–1476. Laufs, H., Holt, J.L., Elfont, R., Krams, M., Paul, J.S., Krakow, K., Kleinschmidt, A., 2006. Where the BOLD signal goes when alpha EEG leaves. NeuroImage 31, 1408–1418. Maes, F., Collingnon, A., Vandermeulen, D., Marchal, G., Suetens, P.,1997. Multimodality image registration by maximization of mutual information. IEEE Tr. MedIm. 16 (2), 187–198. Moosmann, M., Ritter, P., Krastel, I., Brink, A., Thees, S., Blankenburg, F., Taskin, B., Obrig, H., Villringer, A., 2003. Correlates of alpha rhythm in functional magnetic resonance imaging and near infrared spectroscopy. NeuroImage 20, 145–158. Raichle, M.E., Snyder, A.Z., 2007. A default mode of brain function: a brief history of an Evolving idea. NeuroImage 37 (4), 1083–1090. Raichle, M.E., MacLeod, A.M., Snyder, A.Z., Powers, W.J., Gusnard, D.A., Shulman, G.L., 2001. A default mode of brain function. Proc. Natl. Acad. Sci. 98, 676–682. Shmueli, K., van Gelderen, P., de Zwart, J.A., Horovitz, S.G., Fukunaga, M., Jansma, J.M., Duyn, J.H., 2007. Low-frequency fluctuations in the cardiac rate as a source of variance in the resting-state fMRI BOLD signal. Neuroimage 38, 306–320. Tyvaert, L., LeVan, P., Grova, C., Dubleau, F., Gotman, J., 2008. Effects of fluctuating physiological rhythms during prolonged EEG-fMRI studies. Clin. Neurophysiol. 119, 2762–2774.