Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Neuroimage. Author manuscript; available in PMC 2010 June 14.
Published in final edited form as:
PMCID: PMC2884198

Spectral spatiotemporal imaging of cortical oscillations and interactions in the human brain


This paper presents a computationally efficient source estimation algorithm that localizes cortical oscillations and their phase relationships. The proposed method employs wavelet-transformed magnetoencephalography (MEG) data and uses anatomical MRI to constrain the current locations to the cortical mantle. In addition, the locations of the sources can be further confined with the help of functional MRI (fMRI) data. As a result, we obtain spatiotemporal maps of spectral power and phase relationships. As an example, we show how the phase locking value (PLV), that is, the trial-by-trial phase relationship between the stimulus and response, can be imaged on the cortex. We apply the method to spontaneous, evoked, and driven cortical oscillations measured with MEG. We test the method of combining MEG, structural MRI, and fMRI using simulated cortical oscillations along Heschl’s gyrus (HG). We also analyze sustained auditory gamma-band neuromagnetic fields from MEG and fMRI measurements. Our results show that combining the MEG recording with fMRI improves source localization for the non-noise-normalized wavelet power. In contrast, noise-normalized spectral power or PLV localization may not benefit from the fMRI constraint. We show that if the thresholds are not properly chosen, noise-normalized spectral power or PLV estimates may contain false (phantom) sources, independent of the inclusion of the fMRI prior information. The proposed algorithm can be used for evoked MEG/EEG and block-designed or event-related fMRI paradigms, or for spontaneous MEG data sets. Spectral spatiotemporal imaging of cortical oscillations and interactions in the human brain can provide further understanding of large-scale neural activity and communication between different brain regions.

Keywords: Inverse, MEG, Evoked power, Phase locking, Oscillation, Brain, fMRI


Oscillations and synchronization of oscillations (phase-locking) have become fundamental principles in physics and engineering since the first studies by Huygens in 1673. In both magnetoencephalography (MEG) and electroencephalography (EEG), the first signal measured was the alpha rhythm that is thought to be a cortical oscillation originating from the visual cortex and the parietal occipital junction. Due to advances in both source estimation and signal processing techniques, such as wavelet transforms, quantitative analysis of the synchronization of oscillations has recently attracted a lot of interest to understand how information is encoded and transmitted in the brain. A growing number of studies in both humans and animals suggest that cortical oscillations, or rhythmic activity at a variety of frequencies, might be involved in encoding information, and that synchronization (phase-locking) of cortical oscillations might bind and transmit disparate information across brain regions (see recent reviews: (Singer, 1999; Tallon-Baudry and Bertrand, 1999; Varela et al., 2001). Recent biologically realistic models, as well as in vitro and in vivo evidence, suggest that gamma-band oscillations approximately 40 Hz are vital in local information processing, while beta oscillations approximately 20 Hz are used to maintain synchrony over longer distances (Kopell et al., 2000). In the auditory system, cortical oscillations are well documented (Pantev et al., 1991, 1995) and are abnormal in diseased states such as schizophrenia (McCarley et al., 1999). The phase relations in the gamma band have been shown to be important in understanding object perception (Rodriguez et al., 1999; Tallon-Baudry and Bertrand, 1999; Tallon-Baudry et al., 1996, 1997a,b, 1998a,b, 2001), and such relations are abnormal in schizophrenia (Kwon et al., 1999).

Many methods of using spectral analysis on EEG and MEG data have been reported over the past several decades. Some propose wavelet-based techniques to measure neural synchrony between electrodes or between MEG sensors and an external signal, such as an electromyogram (EMG). However, sensors contain information from a variety of brain sources, some of which might overlap. Further, a deep source oriented in the horizontal plane could produce a false apparent synchrony over anterior and posterior sensors. Because source estimates overcome many of these limitations, several methods of source localization of cortical oscillations have been proposed (David et al., 2002; Gross et al., 2001; Hauk et al., 2002; Jensen and Vanni, 2002). Here, we propose a new method that provides the following new features: (1) it creates a map of the cortical oscillations on the cortical surface; (2) it is able to detect rapidly changing power and phase relations; and (3) it can combine prior information from other functional imaging techniques such as fMRI.

The technique presented in this paper is based on a combination of a cortically constrained minimum norm estimate (MNE) (Dale and Sereno, 1993; Dale et al., 2000; Hamalainen and Ilmoniemi, 1984) and wavelet-based spectral analysis employing a complex Morlet wavelet (Tallon-Baudry and Bertrand, 1999; Tallon-Baudry et al., 1997a; Varela et al., 2001). Wavelets preserve a high temporal resolution in the gamma band at approximately 40 Hz necessary to image rapidly time-varying oscillations (Tallon-Baudry and Bertrand, 1999; Tallon-Baudry et al., 1997a; Varela et al., 2001). We evaluate this new technique using both simulated and real MEG data. We validate the technique using a well-characterized simple paradigm: the 40-Hz steady-state-evoked magnetic field. This paradigm was chosen because it produces a well-known periodically varying source in and around Heschl’s gyrus (HG), along the posterior superior temporal gyrus (STG). We can then quantify the accuracy of the localizations of this oscillation using our proposed method. Further, because it is possible to control the magnitude and phase of the stimulus, we can produce a response with a driven phase relationship necessary to evaluate the method. We also demonstrate that this technique reveals spontaneous oscillations, such as the alpha rhythm, in the cortical space.

Materials and methods

Spectral dynamic statistical parametric mapping

Our method to analyze oscillatory activity is called spectral dynamic statistical parametric mapping (sdSPM). The overall method is summarized in Fig. 1. It involves the following computational steps: (1) wavelet transformation MEG sensor data, (2) calculation of the spectral-power minimum-norm estimate (MNE) and the noise-normalized MNE (also called dynamic statistical parametric mapping or dSPM), and (3) calculation of the phase-locking values. To image the induced activity, these steps are repeated for each trial and the results are averaged across trials. If only the spectral power of the evoked response is desired, steps (1) and (2) can be calculated together using the signals averaged over epochs.

Fig. 1
The schematic diagram illustrating the process of using raw MEG data to calculate the phase locking value (PLV) and the frequency specific power, as well as the noise normalized power, of the evoked response on the cortical space with temporal resolution ...

The wavelet transform

To reduce the computational burden, we performed a wavelet transformation of the MEG signals before performing the MNE. The MEG data were first filtered using a continuous wavelet transform (CWT). Temporal convolution of the MEG signal with the wavelet centering at center frequency f and time t produces a temporally localized frequency response


where Ψf,t is the wavelet mother function and the superscript * indicates the complex conjugate. Further, yi(τ) represents the signal in one MEG sensor and wi(t,f) the wavelet coefficients of trial i. Note that for multiple-channel MEG measurements, this convolution is performed over each individual channel yi(τ) separately. In the transformation, we employed the complex Morlet wavelet (Lachaux et al., 1999; Le Van Quyen et al., 2001):


which is the product of a complex exponential of a sinusoidal wavelet, centered at frequency f modulated by a Gaussian with mean τt and variance σ2. To select 40-Hz signals, we set f to be 40 Hz and σ2 to be 7/2 πf.

Source estimation

To localize the current sources underlying the MEG signals, we employed the cortically constrained minimum-norm estimate (MNE) (Dale et al., 2000; Hamalainen et al., 1993). The MNE can be calculated by applying a linear inverse operator W to the measured signals:


where x(t) represents the n-channel MEG data at time t and y(t) are the corresponding current values in the cortical source space. An expression for W is obtained in closed form by minimizing


where C and R denote covariance matrices of the noise and sources, respectively. A is the gain matrix, that is, the solution of the free source orientation forward problem, thus having three columns for each source location, representing the solution according to each orthogonal direction. λ2 is a regularization parameter (Dale et al., 2000; Liu et al., 2002), and ||.||2 indicates the L2 norm. Minimization of Eq. (4) over y yields


We calculate the regularization parameter by


A fixed value of 5 was used for the signal-to-noise ratio (SNR), which reflects the value in many evoked response experiments.

The source covariance matrix incorporates a priori assumption on the spatial distribution of the source currents. It has been proposed that spatial information from fMRI can be incorporated into the MNE by employing a diagonal R with larger elements at the locations of significant fMRI activity (Liu et al., 1998). The MNE is known to have a bias towards superficial currents, caused by the attenuation of the MEG lead fields with increasing source depth. It is possible to compensate for this tendency by scaling the columns of A with a function, denoted here by fk for the kth dipole, increasing monotonically with the source depth. A commonly used choice is


where ap is the pth column of A, the subscript k denotes the location index, and γ is a tunable parameter. Whereas Fuchs et al. (1999) suggest γ = 0.5, we have found in our simulations that this does not provide sufficient compensation. Our preferred value is slightly bigger (γ = 0.6).

The noise covariance matrix C was calculated from the 2-min continuous MEG sensor measurements without the subject (the “empty room noise”) on the same day of the experiment. This is done to prevent biasing the inverse solution against spontaneous brain activation which would inevitably occur with the subject in the measurement device. Here, noise is defined as signals not emanating from the subject.

Spectral power and phase locking value calculation

After the wavelet transformation of the MEG sensor data and subsequent linear transforming by the inverse operator W of each individual trial, we have


We use discrete time index n here to denote the implementation of the calculation in discrete time instants. Here, Ω⃑i (n) denotes the frequency-specific source signal in brain space at time n of the ith trial, and [multiply sign in circle] denotes the one-dimensional convolution on the temporal domain. The convolution is calculated for each channel of yi (n) separately. After source localization, the time-varying power [P with right harpoon above] (n) and synchronization index Θ⃑ (n) were calculated according to



where diag (●) denotes the diagonal entries of the operand in the form of a column vector, and the superscript H denotes the Hermitian, that is, the transpose complex conjugate. Further, Ω⃑R (n) denotes the wavelet coefficient from the convolution of the filtering matrix and the reference time course. Θ⃑ (n) is the phase locking value (PLV) (Lachaux et al., 1999) describing the consistent phase locking between the measurement and the reference time series. In the perfect synchronized case, PLV is 1. PLV is zero when the phase relationship is purely random. In our experiment, Ω⃑R (n) is the 40 Hz sinusoidal oscillation calculated from the complex Morlet wavelet. Note that the power and synchronization index was averaged across trials, not across different time instants.

Noise normalized spectral power using MNE and statistical inference of PLV

Power estimates on the cortical surface were contrasted with the baseline power averages. Note that here the baseline is calculated from the raw data with multiple trials to include the trial-to-trial variability, rather than an averaged baseline as traditionally done with ERP analysis. We adopted the F-statistic, which is the ratio of the instantaneous power at the frequency of interest over the averaged baseline power of the same frequency, to test if statistical significant activities exist (Dale et al., 2000):


where ybaseline(n) denotes the temporal concatenation of baseline measurements from the raw data on MEG sensors, COV (·) denotes the covariance operator, and Cf thus denotes the complex baseline covariance. In the source localization without a priori dipole orientations, under the null hypothesis, [S with right harpoon above]power (n) will follow the F distribution with 3 degrees of freedom for the numerator and a very large number of degrees of freedom for the denominator (the number of trials multiplied by the samples within each trial).

The statistical significance of phase-locking value can also be derived from the Rayleigh test (Babiloni et al., 2002; Fisher, 1993; Simoes et al., 2003). Specifically, when the number of trials in the calculation is large (exceeding 50), the cumulative distribution function of PLV can be approximated by the following:


Fig. 1 summarizes the computation of MNE, noise normalized MNE, and PLV on the cortical surface combing MEG, anatomical MRI, or functional MRI.


We performed a simulation of the sustained 40-Hz activity field emanating from the left superior temporal gyrus (STG), along Heschl’s gyrus (HG). The anatomic information was derived from the structural MRI acquired for the subject of the MEG experiment (see Structural MRI section below for further details). After identifying this region of interest, we placed a patch of current dipoles perpendicular to the cortical surface at HG. The cortical surface geometry was provided from an automatic cortical surface segmentation algorithm (Dale et al., 1999; Fischl et al., 1999). The temporal activity of the source was modeled as a 40-Hz oscillation, with the onset at 50 ms. A sustained 40-Hz activity was maintained for 200 ms. The complete epoch was simulated to be 600 ms with 100 ms baseline. We superimposed Gaussian noise of zero mean and unit variance on each sensor to create an overall signal-to-noise ratio (SNR) for each trial of 0.01 before filtering by the wavelet. Trials (100) were produced to resemble a typical MEG experiment condition. Additionally, the SNRs (before wavelet filtering) were adjusted between 0.1 and 0.001 to test the detection power of different inverse approaches. Prior knowledge of neural activation sites was incorporated by manipulating the R matrix in Eq. (5) for either no fMRI weighting or 90% fMRI weighting. The 90% weighting of the fMRI was chosen to provide the best compromise between separation of activity from correctly localized sources and minimization of error caused by missing sources (Liu et al., 1998). In the case without fMRI weighting, R is an identity matrix, whereas with 90% fMRI weighting, entries corresponding to simulated source were adjusted to be 0.9, and all other entries were adjusted to be 0.1. The MEG signals were calculated with a boundary element model (BEM) (Hamalainen et al., 1993; Oostendorp and van Oosterom, 1989) without spatial decimation, which yields a source space of approximately 340,000 source locations. The inverse operator in the simulation was computed with a spatially decimated forward operator at 7-mm resolution on the cortical surface. The fMRI prior was mapped onto a 7-mm decimated BEM of the gray matter segmented using FreeSurfer ( (Fischl et al., 1999). The location of the simulated source and the simulated fMRI priors is shown in Fig. 2. We used the receiver operating characteristic (ROC) curves to quantify the detection power of the inverse operators. Specifically, true positive rate TPR (α) and false positive rate FPR (α) of the inverse detection × thresholded at level α were calculated as

Fig. 2
The locations of the simulated current sources around Heschl’s gyrus (HG) and the superior temporal gyrus (STG) at the temporal lobe of the left hemisphere are rendered on the inflated cortical surface. The sources are spatially distributed on ...

Here area (Xsim) indicates the area of the simulated sources, area (whole brain) indicates the area of the whole brain, and area (X > α) indicates the area where the inverse solution has value greater than, and ∩ denotes the Boolean AND operator, respectively.


For the MEG and fMRI measurements, a single healthy male volunteer aged 30 was studied. He had no previous history of hearing loss and tested within the normal range with our screening audiogram. All research procedures were performed in accordance with the institutional committee on human research. Written and oral informed consents were both obtained from the volunteer.

MEG acquisition for auditory stimulus

Neuromagnetic responses were recorded using a 306-channel MEG system (Elekta Neuromag OY, Espoo, Finland). During the MEG recording, the subject sat in the chair with his head inside the MEG helmet. Four head position indicator (HPI) coils were placed on the subject’s head to coregister the anatomical MRI and the MEG sensors. The measurement bandwidth was 0.1–172 Hz and the data were digitized at 600 Hz.

A measurement of the room without a subject was made for 4 min that was later used for estimating the noise covariance matrix (see above). After the subject was seated inside the dewar, a 2-min measurement was obtained, while the subject was instructed to rest with eyes open. Another 2-min measurement with the subject’s eyes closed was also made.

Auditory stimuli were delivered by an HP PC driving a Digidesign sound card to Eartone ER3A transducers (Etymotic, Elk Grove, IL, USA). The hearing threshold was defined as the stimulus intensity where the subject responded correctly approximately 70% of the time. Stimuli consisted of a series of 1-ms auditory clicks at 25-ms intervals (40 Hz). Stimuli were presented at 60 dB sensation level (SL). Click trains were generated using Matlab (Mathworks, Natick, MA) and presented independently, with each click burst occurring 90 times. Monaural stimuli were presented alternatively to the left and right ear with 1-s duration and 1 s of silent gap, yielding 4-s trial including 2 s for left ear stimuli (1-s auditory click and 1-s silence) and 2 s for right ear stimuli. The experimental paradigm is shown in Fig. 3.

Fig. 3
The experimental paradigm of auditory stimuli and MEG acquisition. Right and left aural clicks at 40 Hz were sent alternatively every 4 s. Total 90 trials (6 min) were collected for the subject.

MEG acquisition for median nerve stimulus

In the somatosensory study, the right median nerve was stimulated at the wrist with 0.5-ms constant current pulses whose amplitude was just above the motor threshold. The interstimuli interval between current pulses was 4 s. A 306-channel MEG system (Neuromag, Helsinki, Finland) was used to record the neuromagnetic responses. The measurement bandwidth was 0.03–250 Hz and the data were digitized at 1004 Hz. About 100 responses were averaged.

Structural MRI

The brain anatomy was imaged by high-resolution T1-weighted 3D volume MRI using an MPRAGE sequence (TI/TR/TE/flip =1100 ms/2530 ms/3.49 ms/7°, partition thickness = 1.33 mm, matrix = 256 × 256, 128 partitions, field of view = 21 × 21 cm) using a 1.5-T MRI scanner (SIEMENS Medical Solutions, Erlangen, Germany). The anatomical geometry of the gray–white matter surface in the cortex was subsequently derived from an automatic segmentation algorithm to yield a triangulated model with approximate 340,000 vertices (Dale et al., 1999; Fischl et al., 1999). Subsequently, forward matrices were separately calculated using either (1) a decimated a source space of approximately 7000 dipoles with a 7-mm distance between the nearest two dipoles, or (2) using all the vertices in the cortical surface triangulation. The geometry of the inner skull required for the boundary element model was derived from the segmented MRI data with a resolution of approximately 5000 faces (Hamalainen et al., 1993; Oostendorp and van Oosterom, 1989).

Functional for auditory stimulus

The fMRI experimental paradigm is illustrated in Fig. 4. We employed a sparse imaging paradigm (Edmister et al., 1999; Hall et al., 1999) where auditory clicks were presented during the 7-s silent period between consecutive MRI acquisitions to mitigate the acoustic noise generated by gradient switching during the image acquisitions. We used a gradient echo BOLD pulse sequence on a 3-T scanner (SIEMENS Medical Solutions, Erlangen, Germany) with imaging parameters of TR = 9 s, TE = 30 ms, flip angle = 90°, FOV = 200 × 200 mm, matrix = 64 × 64, slice thickness = 5 mm with 10% gap, 24 slices. The analysis of the fMRI data included the motion artifact correction using a 12-parameter rigid-body transformation and spatial smoothing of a 3D Gaussian kernel with full-width-half-maximum (FWHM) of 6 mm. The experimental data were subsequently analyzed by Student’s t test to estimate the statistical significant image voxels between the trials of auditory clicks and the trials of null baseline condition.

Fig. 4
The fMRI experimental paradigm for auditory stimuli. In each trial, seven auditory clicks were present during the 7-s inter-EPI silence period, followed by the 2-s EPI acquisition. Total 50 trials were collected for the subject.



The simulation results are shown in Fig. 5A. We averaged the responses of the 200-ms sustained oscillation period to generate snapshots of the power, noise normalized power, and the PLV on cortical surface. The power estimates were normalized between 0 and 1 to illustrate the relative spatial distribution of the dipole estimates. Noise-normalized power estimates were kept the original scale to reflect the averaged F-statistics (Dale et al., 2000). After 100 trials, all three metrics localized to the correct hemisphere near the temporal lobe region. Without fMRI prior, the power estimate split at the Heschl’s gyrus into two separate activation loci. The 90% fMRI-weighted power estimates provided a well-matched spatial distribution to the location of the simulated cortical patch. The noise-normalized power estimates cover the location of the simulated cortical patch but extend to the insula. However, the 90% fMRI-weighted noise-normalized power estimates generate a phantom activation at the inferolateral aspect of the central sulcus. Phase locking values with and without fMRI priors both cover the location of the simulated patch. Like noise-normalized power estimates, both PLV simulations also show phantom insular activity and even 90% fMRI-weighted PLV shows phantom activation along the inferolateral central sulcus. Fig. 5B shows the receiver operating characteristic (ROC) curves for the spatial sensitivity and specificity of power, noise-normalized power, and PLV inverse using either no fMRI prior or 90% fMRI prior. Note that without fMRI prior, noise normalization can increase sensitivity and specificity of the detection, compared to MNE. Nevertheless, with 90% fMRI weighting, fMRI-weighted MNE outperforms noise normalized power estimates by shifting the ROC curve to the left. Fig. 5C shows the ROC curves when SNR varied between 0.001 and 0.1 (before wavelet filtering). The distribution of ROC curves for PLV and noise-normalized MNE without fMRI prior moved to the left of the distribution of ROC curves for MNE without fMRI prior, indicating that PLV and noise-normalized MNE have higher detection power than MNE when no fMRI priors are available. Provided 90% fMRI prior, fMRI-weighted MNE has the leftmost ROC curve distribution than PLV and noise-normalized MNE over a variety of response SNRs. This validated the advantage of using fMRI prior in MNE to increase the detection power.

Fig. 5Fig. 5
(A) Estimated 40-Hz activity along Heschl’s gyrus in the simulation study. The three rows show the MNE 40-Hz wavelet power (P[var pi], Eq. (9)), the noise-normalized MNE 40-Hz wavelet power (Θ[var pi], Eq. (9)), and the 40-Hz PLV (P ...

Spontaneous MEG measurement and localization

Spontaneous signals were measured in the MEG sensors covering the whole brain including the parietal lobe. Time–frequency analysis using complex Morlet wavelets centered between 6 and 14 Hz revealed that the maximal spectral power detected was approximately 10 Hz. Note that the 10-Hz spontaneous oscillation was not sustained but rather occurred in bursts (see Fig. 6A).

Fig. 6
The most prominent spontaneous responses at 10 Hz during the “eye-closed” condition (A), the auditory stimuli-evoked responses (B), and the somatosensory-evoked responses by median nerve stimulation (C). The top graphs are the signal strength ...

The spontaneous 10-Hz oscillation (alpha) localized bilaterally, predominantly in the right hemisphere. Before noise normalization, the 10-Hz power localized to the superficial parietal gyrus around the occipital parietal junction. Using the noise covariance matrix from 2-min empty room raw data, the noise-normalized 10-Hz power localized the spontaneous activities to retrosplenial cortex at the occipital parietal junction (Fig. 7), consistent with prior imaging studies (Gross et al., 2001; Hari and Salmelin, 1997).

Fig. 7
Medial aspects of the 10-Hz power and the noise normalized 10-Hz power on the inflated cortical surface during the subject’s “eye-closed” condition. These snapshots are temporally averaged across 20 s. Note that the noise normalization ...

Auditory-evoked MEG measurement and localization

Fig. 8 shows the most prominent 6 MEG sensor measurements at the temporal lobe of the left hemisphere. Note that although the sustained response was evident by visual inspection of the waveforms, the time domain signal mixed multiple components including sustained oscillation and primary-evoked response. After wavelet filtering, the frequency analysis revealed a major component of the sustained 40 Hz oscillation, as shown in Fig. 6B. Specifically, at the maximal power channel, we observed the P50 peak followed by N100 peak around 50 and 100 ms after the onset of the auditory stimuli. Subsequently, a sustained oscillation response began approximately 400 ms after stimulus onset and lasted for approximately 400 ms.

Fig. 8
Three x-gradiometer and y-gradiometer MEG sensor measurements around the left hemisphere temporal lobe with the auditory stimulus. In the figure, “grad-x” indicates x-gradiometers, and “grad-y” indicates y-gradiometers. ...

Using auditory clicks, fMRI revealed that BOLD activation along the left Heschl’s gyrus and the STG contrasted with the silent baseline. This is illustrated in Fig. 9A at a t statistics threshold of 12 (uncorrected P < 10−10). The decimated dipoles (with 7-mm spacing) coincided with the fMRI active region were recruited as the priors for MNE. The source localization of 40-Hz power, noise-normalized 40-Hz power, and phase-locking values with 90% or without fMRI priors were shown in Fig. 8B. In the power illustrations, we average the power between 500 and 900 ms (see Fig. 9B), and subsequently normalized between 0 and 1 to illustrate the brain areas with power exceeding 30% of the maximum. Similar temporal averaging and source estimate scaling were done for noise-normalized power and PLV calculation. Note that without fMRI weighting, the estimation of source power shifted toward the STG. Using fMRI weighting of 90%, significant power was centered in and around Heschl’s gyrus as indicated by fMRI priors. For PLV localizations, inverse estimates with or without fMRI priors do not differ significantly. Comparing power and PLV localizations, we found that PLVs with or without fMRI priors can detect stimulus-sensitive regions including Heschl’s gyrus.

Fig. 9
(A) The t statistics from three runs of fMRI experiments with paradigm shown in Fig. 2. The t statistics were calculated between the auditory click condition and the silence condition. The areas with negative l0-based logarithm of the P values (uncorrected ...

Somatosensory-evoked MEG measurement and localization

Fig. 6C shows the time domain and spectral domain somatosensory-evoked response by the median nerve stimulation at the MEG sensor close to the central sulcus of right hemisphere. In the time–frequency representation (bottom of Fig. 6C), we observed a transient response concentrating at 17 Hz, consistent with the upper portion of the mu-rhythm (Hari and Salmelin, 1997). The transient response begins 40 ms after median nerve stimulation, and lasts for approximately 70 ms. Using wavelet centered at 17 Hz and performing the MNE, we can localize this transient oscillation to the right hemisphere central sulcus. This is shown in Figs. 10A and B, which are MNE and noise-normalized MNE 17-Hz power, respectively. Using a dipole centered in the right primary sensory cortex (SI) with maximal power estimated by MNE as the reference, we calculated the PLV and observed the synchronous activity at contralateral hemisphere (Fig. 10C). Significant phase locking phenomena were observed at the left hemisphere central sulcus SI area and left hemisphere SII area, consistent with other MEG reports (Simoes et al., 2003).

Fig. 10
The 17-Hz wavelet-based spectral spatiotemporal maps using MNE and noise-normalized MNE at the right hemisphere (left-hand median nerve stimulation). Based on the maximal power dipole at right hemisphere around central sulcus, a left-hemisphere PLV map ...


Here, we demonstrated a computationally efficient method of a wavelet-based, anatomically constrained MNE to create an image of signal frequency power and phase on the cortical surface. In particular, we used individual anatomy (head and brain) to further guide the source estimation, and display it on the cortex. The technique visualizes short-lasting activity including both the spectral power over some frequency band and the phase relationship, either between brain regions or between a brain region and the stimulus.

We emphasize the importance of choosing and calculating the appropriate noise covariance matrix. For spontaneous activities, such as 10-Hz alpha rhythm in our results, the noise refers to the noise observed extraneous from the subject’s brain activity. Therefore, empty room recordings were used to calculate this noise covariance matrix.

For the noise normalization, the term “noise covariance” actually refers to the covariance of the prestimulus interval, that is, the baseline covariance. An incorrect choice of the baseline covariance causes either overestimation or underestimation of the baseline activity as well as the possibility of mislocalization of the statistical significance calculated by noise normalization. In addition, the baseline covariance should always be computed from the raw data to increase the number of samples available to the estimator. A baseline covariance matrix estimated from averaged data always undercharacterizes the variability of the measurement and thus leads to overestimated significance in the statistical tests. This is true because the covariance of the averaged data does not account for trial-to-trial variations. Therefore, raw data in appropriate “baseline” conditions have to be chosen appropriately for noise normalization.

Comparing our simulations and the results from real data, we validated that the reported technique can be used to detect frequency-specific responses at signal to noise ratio (SNR) as small as 0.01, if up to 100 trials of measurements are available and the activity is phase-locked to the stimulus (evoked). However, if the activity is not phase locked to the stimulus (induced), the SNR must be greater than 1, a limitation that is true for measuring induced spectral power. Nevertheless, the magnitude of the F statistic and PLV differed in the simulations and the actual MEG measurements. Potentially, this may be due to (1) anatomical variations in the neural sources (not precisely perpendicular orientations of dipole sources to the local cortical surface), (2) anatomical information errors from MRI, and (3) time-varying activities in real brain responses (changes of spatial distribution and temporal duration of the oscillations). In the Bayesian estimation theory framework, this can be further improved by providing more accurate prior information regarding both the spatial and temporal properties of the brain activations, such as more accurate quantitative description between BOLD fMRI and MEG measurements. This is a current research topic our laboratory is actively investigating. It is hoped that this method will make such research possible.

Comparison with other source estimation methods

Currently, several methods have been proposed to map cortical oscillations, each of which has a particular strength for its application, but none of which meets all the requirements proposed here. Minimum current estimate (MCE), for example (Matsuura and Okabe, 1995; Uutela et al., 1999), assumes sources with exponential probability distribution function. Specifically, Jensen and Vanni (2002) proposed to use MCE to measure simultaneous cortical oscillations, based on the minimum L1-norm estimate, using complex Fourier components from sections of raw MEG signals. Although it is limited to the analysis of spontaneous cortical oscillations, the method could be extended to include measurements of synchronization based on the phase of the Fourier components. It is, however, limited in that it does not have a high-temporal resolution to measure fast changing oscillations and requires a high SNR. Source localizations using beamformer techniques (Van Veen et al., 1997; Vrba et al., 1999) have been also applied to study neural synchrony, such as dynamic imaging of coherent sources (DICS) (Gross et al., 2001). This innovative technique is a promising method of detecting synchronizations between neural networks in the brain, based on source estimation using a spatial filter of coherence. Like the MCE technique, it is limited, however, in that it requires a relatively high SNR and makes the assumption of stationary data, and consequently does not measure short lasting synchronizations. Hashimoto et al. (2001), Sekihara et al. (2001), and Sekihara et al. (2002) recently demonstrated that the stationary data assumption is not necessarily strict.

The use of the MNE to image cortical oscillations has also been proposed by others (David et al., 2002; Hauk et al., 2002). Hauk et al. (2002) compared the MNE to a current source density method using EEG data sets with a wavelet-based frequency domain analysis. Others (David et al., 2002; Lachaux et al., 2002) have recently performed simulations on various versions of the MNE. They concluded that MNE is a preferred method compared to these other minimum L2-norm inverse approaches. Our method differs in several ways. First, our solution computes single trial epochs, or uses clips of spontaneous activity, to create an averaged power spatiotemporal map. Second, we use a BEM to use the features of the actual head of the individual subject’s MRI to create a power sources estimate directly on the cortical surface, utilizing a cortical constraint. Third, we also include the ability to create a phase synchronization map of the brain on the cortical surface. In this case, we show how the phase relation of the stimulus to the response can be created, but the same principle can be used to create maps showing synchronization between brain regions, provided the reference temporal course from one specific brain region. And last, the utilization of prior neural activation information from fMRI can help for correct power localization. This is especially important during the localization of the evoked auditory responses because of highly convoluted cortical surface at the temporal lobe. The anatomy challenges the traditional MNE to give focal and precise source estimations. Without fMRI priors, it is easy to lead to superficial source estimates away from Heschl’s gyrus.

Our method (sdSPM) has similarities to that proposed by David et al. (2002) in that both source estimation methods use a wavelet-based spectral analysis and a MNE. Our method differs in that we calculate the wavelet-based spectral power based on the power in source space (on the brain surface) across trains, instead of calculating the single trial power and single trial phase locking maps as used by David et al. We have also proposed a method of producing a spectral dSPM with noise normalization of the MNE.


Several limitations of this technique are illustrated in this report. The algorithm’s localizing ability is often limited, especially with activity in and around the Sylvian fissure, a common problem with imaging peri-Sylvian activity generally and auditory activity specifically. We can see in our results that probable ‘phantom’ sources are seen in the posterior inferior frontal lobe, insula, and the middle temporal gyrus. These sources, even with actual measurements, are likely phantom as they are seen in the simulations as well due to the close anatomical proximity to Heschl’s gyrus. Due to the anatomical convolution of cortical surface, traditional MNE without a priori activation loci may lead to mislocalization of auditory neural activities. Using fMRI, which provides a higher spatial resolution brain activation estimate, it is possible to constrain the localization process to avoid overestimation of the power of neuromagnetic sources in traditional MNE. For phase-locking value and noise-normalized power estimates, it appears that a priori information could be less essential. This is because of the intrinsic normalization of PLV calculation by disregarding the amplitude of dipole estimates in an individual trial. The noise normalization also compensates for the biases on the magnitude of the estimates. Therefore, PLV and dSPM localization shows similar sensitivity for both superficial sources on the gyri and deep sources in the sulci. Another limitation of our approach is that no deep gray nuclei activity is included in the source model. Thalamocortical interaction is probably responsible for the generation and maintenance of oscillations in the neocortex. Our group is actively investigating this issue and it will be a focus for an upcoming study.


Although both PLV and noise-normalized power can correct biases toward superficial cortical activities, potential overcorrection may occur, as we found in both the simulations and real data that the insula was probably falsely activated. For simultaneous high sensitivity and specificity, fMRI-weighted MNE can be a candidate tool for auditory response localization.

Future applications

The new method detailed here is a general framework and can be extended to include other imaging technologies including PET, optical imaging, and transcranial Doppler. Any functional information that can be localized can be included in the R matrix and used as prior information and thus improve the localizing ability of the MEG or EEG data.


This work was supported by National Institutes of Health Grants R01 HD040712, R01 NS037462, P41 RR14075, the Mental Illness and Neuroscience Discovery Institute (MIND), and the American Heart Association.


  • Babiloni C, Babiloni F, Carducci F, Cincotti F, Rosciarelli F, Arendt-Nielsen L, Chen AC, Rossini PM. Human brain oscillatory activity phase-locked to painful electrical stimulations: a multi-channel EEG study. Hum Brain Mapp. 2002;15:112–123. [PubMed]
  • Dale AM, Sereno MI. Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach. J Cogn Neurosci. 1993;5:162–176. [PubMed]
  • Dale AM, Fischl B, Sereno MI. Cortical surface-based analysis: I. Segmentation and surface reconstruction. Neuroimage. 1999;9:179–194. [PubMed]
  • Dale AM, Liu AK, Fischl BR, Buckner RL, Belliveau JW, Lewine JD, Halgren E. Dynamic statistical parametric mapping: combining fMRI and MEG for high-resolution imaging of cortical activity. Neuron. 2000;26:55–67. [PubMed]
  • David O, Garnero L, Cosmelli D, Varela FJ. Estimation of neural dynamics from MEG/EEG cortical current density maps: application to the reconstruction of large-scale cortical synchrony. IEEE Trans Biomed Eng. 2002;49:975–987. [PubMed]
  • Edmister WB, Talavage TM, Ledden PJ, Weisskoff RM. Improved auditory cortex imaging using clustered volume acquisitions. Hum Brain Mapp. 1999;7:89–97. [PubMed]
  • Fischl B, Sereno MI, Dale AM. Cortical surface-based analysis: II. Inflation, flattening, and a surface-based coordinate system. Neuroimage. 1999;9:195–207. [PubMed]
  • Fisher NI. Statistical Analysis of Circular Data. Cambridge University Press; New York: 1993.
  • Fuchs M, Wagner M, Kohler T, Wischmann HA. Linear and nonlinear current density reconstructions. J Clin Neurophysiol. 1999;16:267–295. [PubMed]
  • Gross J, Kujala J, Hamalainen M, Timmermann L, Schnitzler A, Salmelin R. Dynamic imaging of coherent sources: studying neural interactions in the human brain. Proc Natl Acad Sci U S A. 2001;98:694–699. [PubMed]
  • Hall DA, Haggard MP, Akeroyd MA, Palmer AR, Summerfield AQ, Elliott MR, Gurney EM, Bowtell RW. “Sparse” temporal sampling in auditory fMRI. Hum Brain Mapp. 1999;7:213–223. [PubMed]
  • Hamalainen M, Ilmoniemi R. Helsinki University of Technology, Technical Report TKK-F-A559. 1984. Interpreting measured magnetic fields of the brain: estimates of current distributions.
  • Hamalainen M, Hari R, Ilmoniemi R, Knuutila J, Lounasmaa O. Magnetoencephalography—Theory, instrumentation, and application to noninvasive studies of the working human brain. Rev Mod Phys. 1993;65:413–497.
  • Hari R, Salmelin R. Human cortical oscillations: a neuromagnetic view through the skull. Trends Neurosci. 1997;20:44–49. [PubMed]
  • Hashimoto I, Kimura T, Iguchi Y, Takino R, Sekihara K. Dynamic activation of distinct cytoarchitectonic areas of the human SI cortex after median nerve stimulation. Neuroreport. 2001;12:1891–1897. [PubMed]
  • Hauk O, Keil A, Elbert T, Muller MM. Comparison of data transformation procedures to enhance topographical accuracy in time-series analysis of the human EEG. J Neurosci Methods. 2002;113:111–122. [PubMed]
  • Jensen O, Vanni S. A new method to identify multiple sources of oscillatory activity from magnetoencephalographic data. Neuroimage. 2002;15:568–574. [PubMed]
  • Kopell N, Ermentrout GB, Whittington MA, Traub RD. Gamma rhythms and beta rhythms have different synchronization properties. Proc Natl Acad Sci U S A. 2000;97:1867–1872. [PubMed]
  • Kwon JS, O’Donnell BF, Wallenstein GV, Greene RW, Hirayasu Y, Nestor PG, Hasselmo ME, Potts GF, Shenton ME, McCarley RW. Gamma frequency-range abnormalities to auditory stimulation in schizophrenia. Arch Gen Psychiatry. 1999;56:1001–1005. [PMC free article] [PubMed]
  • Lachaux JP, Rodriguez E, Martinerie J, Varela FJ. Measuring phase synchrony in brain signals. Hum Brain Mapp. 1999;8:194–208. [PubMed]
  • Lachaux JP, Lutz A, Rudrauf D, Cosmelli D, Le Van Quyen M, Martinerie J, Varela F. Estimating the time-course of coherence between single-trial brain signals: an introduction to wavelet coherence. Neurophysiol Clin. 2002;32:157–174. [PubMed]
  • Le Van Quyen M, Foucher J, Lachaux J, Rodriguez E, Lutz A, Martinerie J, Varela FJ. Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony. J Neurosci Methods. 2001;111:83–98. [PubMed]
  • Liu AK, Belliveau JW, Dale AM. Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations. Proc Natl Acad Sci U S A. 1998;95:8945–8950. [PubMed]
  • Liu AK, Dale AM, Belliveau JW. Monte Carlo simulation studies of EEG and MEG localization accuracy. Hum Brain Mapp. 2002;16:47–62. [PubMed]
  • Matsuura K, Okabe Y. Selective minimum-norm solution of the biomagnetic inverse problem. IEEE Trans Biomed Eng. 1995;42:608–615. [PubMed]
  • McCarley RW, Wible CG, Frumin M, Hirayasu Y, Levitt JJ, Fischer IA, Shenton ME. MRI anatomy of schizophrenia. Biol Psychiatry. 1999;45:1099–1119. [PMC free article] [PubMed]
  • Oostendorp TF, van Oosterom A. Source parameter estimation in inhomogeneous volume conductors of arbitrary shape. IEEE Trans Biomed Eng. 1989;36:382–391. [PubMed]
  • Pantev C, Makeig S, Hoke M, Galambos R, Hampson S, Gallen C. Human auditory evoked gamma-band magnetic fields. Proc Natl Acad Sci U S A. 1991;88:8996–9000. [PubMed]
  • Pantev C, Bertrand O, Eulitz C, Verkindt C, Hampson S, Schuierer G, Elbert T. Specific tonotopic organizations of different areas of the human auditory cortex revealed by simultaneous magnetic and electric recordings. Electroencehalogr Clin Neurophysiol. 1995;94:26–40. [PubMed]
  • Rodriguez E, George N, Lachaux JP, Martinerie J, Renault B, Varela FJ. Perception’s shadow: long-distance synchronization of human brain activity. Nature. 1999;397:430–433. [PubMed]
  • Sekihara K, Nagarajan SS, Poeppel D, Marantz A, Miyashita Y. Reconstructing spatio-temporal activities of neural sources using an MEG vector beamformer technique. IEEE Trans Biomed Eng. 2001;48:760–771. [PubMed]
  • Sekihara K, Nagarajan SS, Poeppel D, Marantz A. Performance of an MEG adaptive-beamformer technique in the presence of correlated neural activities: effects on signal intensity and time-course estimates. IEEE Trans Biomed Eng. 2002;49:1534–1546. [PubMed]
  • Simoes C, Jensen O, Parkkonen L, Hari R. Phase locking between human primary and secondary somatosensory cortices. Proc Natl Acad Sci U S A. 2003;100:2691–2694. [PubMed]
  • Singer W. Neuron. Vol. 24. 1999. Neuronal synchrony: a versatile code for the definition of relations? pp. 49–65.pp. 111–125. [PubMed]
  • Tallon-Baudry C, Bertrand O. Oscillatory gamma activity in humans and its role in object representation. Trends Cogn Sci. 1999;3:151–162. [PubMed]
  • Tallon-Baudry C, Bertrand O, Delpuech C, Pernier J. Stimulus specificity of phase-locked and non-phase-locked 40 Hz visual responses in human. J Neurosci. 1996;16:4240–4249. [PubMed]
  • Tallon-Baudry C, Bertrand O, Delpuech C, Permier J. Oscillatory gamma-band (30–70 Hz) activity induced by a visual search task in humans. J Neurosci. 1997a;17:722–734. [PubMed]
  • Tallon-Baudry C, Bertrand O, Wienbruch C, Ross B, Pantev C. Combined EEG and MEG recordings of visual 40 Hz responses to illusory triangles in human. Neuroreport. 1997b;8:1103–1107. [PubMed]
  • Tallon-Baudry C, Bertrand O, Peronnet F, Pernier J. Induced gamma-band activity during the delay of a visual short-term memory task in humans. J Neurosci. 1998;18:4244–4254. [PubMed]
  • Tallon-Baudry C, Bertrand O, Fischer C. Oscillatory synchrony between human extrastriate areas during visual short-term memory maintenance. J Neurosci. 2001;21:RC177. [PubMed]
  • Uutela K, Hamalainen M, Somersalo E. Visualization of magnetoencephalographic data using minimum current estimates. Neuroimage. 1999;10:173–180. [PubMed]
  • Van Veen BD, van Drongelen W, Yuchtman M, Suzuki A. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans Biomed Eng. 1997;44:867–880. [PubMed]
  • Varela F, Lachaux JP, Rodriguez E, Martinerie J. The brainweb: phase synchronization and large-scale integration. Nat Rev, Neurosci. 2001;2:229–239. [PubMed]
  • Vrba J, Anderson G, Betts K, Burbank M, Cheung T, Cheyne D, Fife A, Govorkov S, Habib F, Haid G, Haid V, Hoang T, Hunter C, Kubik P, Lee S, McCubbin J, McKay J, McKenzie D, Nonis D, Paz J, Reihl E, Ressel D, Robinson S, Schroyen C, Sekatchev I, Spear P, Taylor B, Tillotson M, Sutherlilng W. 151-channel whole cortex MEG system for seated or supine positions. 11th International Conference on Biomagnetism.1999.