www.elsevier.com/locate/ynimg
NeuroImage 40 (2008) 1686 – 1700

Five-dimensional neuroimaging: Localization of the time–frequency
dynamics of cortical activity
Sarang S. Dalal,a,b Adrian G. Guggisberg,a,g Erik Edwards,a,e Kensuke Sekihara,fAnne M. Findlay,a Ryan T. Canolty,e Mitchel S. Berger,c Robert T. Knight,eNicholas M. Barbaro,c,d Heidi E. Kirsch,d and Srikantan S. Nagarajan a,b,⁎aBiomagnetic Imaging Laboratory, Department of Radiology, University of California, San Francisco, CA 94143-0628, USA
Joint Graduate Group in Bioengineering, University of California, San Francisco and University of California, Berkeley, USAc
Department of Neurological Surgery, University of California, San Francisco, CA 94143-0112, USAd
Department of Neurology, University of California, San Francisco, CA, 94143-0138, USAeHelen Wills Neuroscience Institute and Department of Psychology, University of California, Berkeley, CA 94720-3190, USAf
Department of Systems Design and Engineering, Tokyo Metropolitan University, Tokyo 191-0065, Japang
Department of Neurology, University of Bern, Inselspital, 3010 Bern, Switzerlandb

Received 28 March 2007; revised 8 January 2008; accepted 17 January 2008
Available online 31 January 2008

The spatiotemporal dynamics of cortical oscillations across humanbrain regions remain poorly understood because of a lack ofadequately validated methods for reconstructing such activity fromnoninvasive electrophysiological data. In this paper, we present a novel
adaptive spatial filtering algorithm optimized for robust source time–
frequency reconstruction from magnetoencephalography (MEG) andelectroencephalography (EEG) data. The efficacy of the method isdemonstrated with simulated sources and is also applied to real MEG
data from a self-paced finger movement task. The algorithm reliably
reveals modulations both in the beta band (12–30 Hz) and high gamma
band (65–90 Hz) in sensorimotor cortex. The performance is validated
by both across-subjects statistical comparisons and by intracranialelectrocorticography (ECoG) data from two epilepsy patients. Interestingly, we also reliably observed high frequency activity (30–300 Hz)
in the cerebellum, although with variable locations and frequencies
across subjects. The proposed algorithm is highly parallelizable andruns efficiently on modern high-performance computing clusters. Thismethod enables the ultimate promise of MEG and EEG for fivedimensional imaging of space, time, and frequency activity in the brain
and renders it applicable for widespread studies of human cortical
dynamics during cognition.
© 2008 Elsevier Inc. All rights reserved.

⁎ Corresponding author. Biomagnetic Imaging Laboratory, Department ofRadiology, University of California, San Francisco, CA 94143-0628, USA.
E-mail address: sri@radiology.ucsf.edu (S.S. Nagarajan).
Available online on ScienceDirect (www.sciencedirect.com).
1053-8119/$ - see front matter © 2008 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2008.01.023

Introduction
Magnetoencephalography (MEG) and electroencephalography(EEG) are functional neuroimaging techniques with millisecond
time resolution (Hämäläinen et al., 1993). Traditionally, MEG andEEG have been used to study evoked responses, i.e., activity that isboth time-locked and phase-locked to a stimulus or task. These
analyses assume a model of neural activity in which responses are
additive and/or phases are reset (Hanslmayr et al., 2007). However,it has been well-known that ongoing MEG/EEG oscillations can besuppressed in response to a stimulus or task since the earliest EEG
research (Berger, 1930); this possibility is not accounted for by theevoked model. Furthermore, the across-trial jitter inherent inresponses to even simple stimuli have been shown to be sufficient
to markedly reduce the amplitude of averaged responses (Michalewski et al., 1986); this effect becomes even more pronounced forhigher frequency bands. Averaging also assumes trial-to-trial phaselocking, which may not be valid for many complex cognitive
paradigms.
Another approach to interpreting MEG and EEG data is toquantify oscillatory aspects of the signals using time–frequency
methods. Typically, modulations of oscillatory activity are describedas event-related spectral power changes (Pfurtscheller and Aranibar,
1977; Pfurtscheller and Neuper, 1992; Makeig, 1993). By comparing the power of neural activity to a quiescent baseline, these types ofanalyses reveal induced responses, i.e., activity that is time-locked
but not necessarily phase-locked. Additionally, the power change
may be negative, termed an event-related desynchronization (ERD),
or positive, termed an event-related synchronization (ERS).
Analyses of ERD and ERS overcome many of the limitations of
evoked response analyses. However, most MEG/EEG time–

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

frequency analyses are conducted on the sensor signals without
source localization, providing only vague information as to whichbrain structures generated the activity of interest.
Several source reconstruction algorithms, each employing a
different set of assumptions, have been proposed to overcome theill-posed inverse problem. Source reconstructions from MEG data can
be classified as either parametric or tomographic. Parametric methods
include equivalent current dipole (ECD) fitting techniques; they often
require knowledge about the number of sources and their approximate
locations and poorly modeled sources with a large spatial extent.
Tomographic methods reconstruct source activity at each voxel(3-D location) in the brain. Spatial filtering techniques avoid thehigh number of parameters and the nonlinear iterative searchrequired by ECD analysis. Nonadaptive spatial filtering techniques,
which include minimum-norm-based methods such as sLORETA(Pascual-Marqui, 2002), use sensor geometry to construct theweights for the spatial filter. Adaptive techniques, on the other hand,additionally use sensor data to create a custom filter depending on
signal characteristics. It has been shown that a class of adaptive
spatial filters known as beamformers (Van Veen and Buckley, 1988)
have the best spatial resolution and performance amongst existing
tomographic methods (Darvas et al., 2004; Sekihara et al., 2005).
Spatial filtering methods have the potential to compute electromagnetic source images in both the time and frequency domains(Robinson and Vrba, 1999; Gross et al., 2001; Sekihara et al., 2001;
Jensen and Vanni, 2002; Dalal et al., 2004). Techniques such as thesynthetic aperture magnetometry (SAM) beamformer have been
employed to examine either the time course of neural sources or the
spatial distribution of power within a specific frequency band
(Robinson and Vrba, 1999). However, published studies typically
employ SAM to generate static fMRI-style images using a large
bandwidth and wide time window—effectively discarding the
temporal resolution advantage of magnetoencephalography. Only afew studies have attempted time–frequency analysis in source space(Singh et al., 2002; Cheyne et al., 2003; Brookes et al., 2004; Gaetz
and Cheyne, 2006; Jurkiewicz et al., 2006). These reports describe amethod in which a single set of beamformer weights are firstcomputed over a wide time window and frequency range; time–
frequency decompositions are then computed from the reconstructed
time series for a few locations of interest. However, as we show in thispaper, weights computed from unfiltered or wideband data may be
inherently biased towards resolving low-frequency brain activity due
to the power law of typical electrophysiological data. Additionally,
responses of shorter duration or outside the fixed time window used
to generate the weights may not be adequately captured.
In this paper, we propose a novel adaptive spatial filtering algorithm
that is optimized for time–frequency source reconstructions from
MEG/EEG data. Performance of this algorithm will first be evaluatedwith simulated data. Then we will demonstrate the method with real
finger movement data, validated with group statistics and intracranialrecordings. The proposed algorithm enables accurate reconstruction of
five-dimensional brain activity from MEG and EEG data, thereby
realizing the ultimate promise of MEG- and EEG-based neuroimaging.
Methods
Definitions and problem formulation
Throughout this paper, plain italics indicate scalars, lowercase
boldface italics indicate vectors, and uppercase boldface italics
indicate matrices.

1687

We define the magnetic field measured by the mth detector coilat time t as bm(t) and a column vector b(t) ≡ [b1(t), b2(t), …, bM(t)]T
as a set of measured data, where M is the total number of detector
coils and the superscript T indicates the matrix transpose. The
second-order moment matrix of the measurement is denoted R, i.e.,Ruhbðt ÞbT ðt Þi, where hi indicates the ensemble average overtrials. When hbðt Þi ¼ 0, R is also equal to the sample covariancematrix. In practice, the covariance is estimated over a subset oflatencies, t ≡ [t1, t2,…, tN], that represents samples from a desired
time window of length N. Defining B(t) ≡ [b(t1), b(t2),…, b(tN)], thecovariance estimate then becomes Rðt ÞuhBðt ÞBT ðt Þi.
We assume that the sensor data arises from elemental dipoles
at each spatial location r, represented by a 3-D vector such that r =(rx, ry , rz). The orientation of each source is defined as a vector
η(r) ≡ [βx, βy, βz], where βx, βy, and βz are the angles between the
moment vector of the source and the x, y, and z axes, respectively.
We define lζm(r) as the output of the mth sensor that would be
induced by a unit–magnitude source located at r and pointing in theζ direction. The column vector lζ(r) is defined as lζ(r) ≡ [lζ1(r),lζ2(r), …, lζM(r)]T. The lead field matrix, which represents thesensitivity of the whole sensor array at r, is defined as L(r) ≡ [lx(r),ly(r), lz(r)]. The lead field vector for a unit-dipole oriented in the
direction η is defined as l(r, η) where l(r, η) ≡ L(r)η(r).
Conventional adaptive spatial filtering
This section reviews an adaptive spatial filter called theminimum variance (MV) scalar beamformer, also referred to as
the synthetic aperture magnetometry (SAM) beamformer (Robinson and Vrba, 1999). An adaptive spatial filter estimate of thesource moment ŝ (r, t) is given byŝðr; t Þ ¼ wT ðrÞbðt Þ

ð1Þ

where w(r) is the weight vector.
The MV scalar beamformer weight vector w(r) is calculated byminimizing wT(r)R(t)w(r) subject to lT(r, η)w(r) = 1. The solution is
known to be (Robinson and Vrba, 1999):R1 ðt Þl ðr; ηÞw ðr Þ ¼ T:
l ðr; ηÞR1 ðt Þl ðr; ηÞð2Þ

Finally, in the absence of a priori orientation information from,e.g., MRI, an optimal orientation ηopt(r) must be determined. The
typical approach to determining ηopt is to compute the solution thatmaximizes output power with respect to η (Sekihara and Scholz,
1996). Our approach is to compute the solution that maximizes
output SNR (Sekihara et al., 2004):l T ðr; ηÞR1 ðt Þl ðr; ηÞηopt ðrÞ ¼max Tηl ðr; ηÞR2 ðt Þl ðr; ηÞð3Þ

As shown in the study of Sekihara et al. (2004), the solutionfor ηopt is v3, the eigenvector corresponding to the minimum
eigenvalue of: T1  T

L ðrÞR1 ðt ÞLðrÞ
L ðrÞR2 ðt ÞLðrÞ vj ¼ gj vj ;
ð4ÞThe estimated source power P̂s(r, t) can be computed from theweights w and covariance R(t):P̂s ðr; t Þhŝðr; t Þ2 i ¼ h wT ðrÞBðt Þ BT ðt ÞwðrÞ i ¼ wT ðrÞRðt ÞwðrÞ
ð5Þ

1688

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

The sensor noise power σ2(t) may be obtained from calibration
measurements of the MEG system or estimated by computing theminimum eigenvalue of R(t). Then, the power of projected sensor
noise P̂N may be estimated by replacing R(t) with σ2(t)I:P̂N ðrÞ ¼ wT ðrÞ r2 ðt ÞI wðrÞ ¼ r2 ðt ÞwT ðrÞwðrÞ
ð6Þ
Often, one is interested in the change in power from a control(i.e., baseline) time window to an active time window, i.e., a dualcondition paradigm. These windows are denoted as vectors of timesamples, tcon and tact, respectively. In this case:P̂con ðrÞ ¼ P̂s ðr; t con Þ ¼ wT ðrÞRcon wðrÞð7ÞP̂act ðrÞ ¼ P̂s ðr; t act Þ ¼ wT ðrÞRact wðrÞð8ÞFinally,
FdB ðr; n; f Þ ¼ 10 log10P̂act ðr; n; f Þ  P̂N ðr; n; f Þ:P̂con ðr; n; f Þ  P̂N ðr; n; f Þ

ð13Þ

However, while spectrograms may be constructed from theconventional beamformer in this fashion, the weights are still optimized for the wide tact and tcon windows used to compute w(r).
MEG/EEG spectra follow the power law, implying that weights
generated from unfiltered data are inherently biased towards lowfrequency activity.
Frequency-dependent weight computation

where Rcon ≡ R(tcon), the covariance of the control window, andRact ≡ R(tact), the covariance of the active window.
In order to improve numerical stability and ensure an appropriately matched baseline period, the same orientation ηopt(r) andw(r) must be used to compute P̂act(r) and P̂con(r). This ensures thatthe magnitude of sources are comparable between the active andcontrol periods; it also decreases the likelihood of resolving false
sources. Thus, ηopt(r) and w(r) may be computed using the averagecovariance of the active and control periods, i.e., by substituting R =(Ract + Rcon)/2. Note that tcon must be the same length as tact.
The contrast between P̂act and P̂con can then be expressed as a
pseudo-t difference P̂act − P̂con or an F-ratio P̂act / P̂con. If thecontribution of projected sensor noise is subtracted, the ratio
becomes F = (P̂act − P̂N) / (P̂con − P̂N). In this paper, we will use thenoise-corrected F-ratio expressed in units of decibels:Pˆ  Pˆ N:
FdB ¼ 10 log10 ˆ act
Pcon  Pˆ N

ð9Þ

Time–frequency extension of conventional beamformers
It is often desirable to compute contrasts for multiple activation
windows and possibly multiple baseline windows, relative to
specific experimental or cognitive events. The resulting contrasted
spectrogram is a time–frequency representation of source events. In
order to obtain such a representation from the conventional
beamformer, one may directly compute the spectrogram of thesource time series from Eq. (1), contrasting it with the spectrogramof the control period (Singh et al., 2002; Cheyne et al., 2003).
Another approach is to apply the weights w(r) computed above—with R estimated from long time windows tact and tcon spanning theentire duration of interest—to a new set of covariance estimates
generated from filtered and segmented data. First, the data is passed
through a filter bank and partitioned into several overlapping active segments, τact[n], and a control segment, τcon, where the subscript n refers to the index of the time window. (These windows
are shorter than tact and tcon.) Then, covariances are computed for eachresulting time–frequency window, yielding R̃ act(n, f ) ≡ R(τact[n], f )and R̃con( f ) ≡ R(τcon, f ), where f corresponds to the index of thefrequency band. Power maps may be computed directly by replacing
Ract and Rcon with R̃ act(n, f ) and R̃ con( f ), respectively:P̂con ðr; n; f Þ ¼ wT ðrÞ R̃con ð f ÞwðrÞ

ð10Þ

P̂act ðr; n; f Þ ¼ wT ðrÞ R̃act ðn; f ÞwðrÞ

ð11Þ

P̂N ðr; n; f Þ ¼ r2 ðn; f ÞwT ðrÞwðrÞ

ð12Þ

Therefore, in order to better resolve low-amplitude, highfrequency activity, one approach is to calculate a different set ofweights for each frequency band:R1 ð f Þl ðr; ηÞwðr; f Þ ¼ T
l ðr; ηÞR1 ð f Þl ðr; ηÞð14Þ

where R( f ) is the sample covariance matrix generated from B(t)filtered for the frequency band of interest, and η = ηopt(r, f ), i.e.,the optimum orientation computed using R( f ). The corresponding
power at each voxel for each frequency band is:P̂s ðr; f Þ ¼ wT ðr; f ÞRð f Þwðr; f Þ

ð15Þ

Again, the powers of an active window and a control window
may be computed as follows:P̂con ðr; f Þ ¼ P̂s ðr; t con ; f Þ ¼ wT ðr; f ÞRcon ð f Þwðr; f Þð16Þ

P̂act ðr; f Þ ¼ P̂s ðr; t act ; f Þ ¼ wT ðr; f ÞRact ð f Þwðr; f Þ

ð17Þ

The time–frequency representation may be computed either
from the source time series, or, as shown here, by using w(r, f )from Eq. (14) and replacing R( f ) from Eq. (15) with R̃act(n, f ) andR̃con( f ):P̂con ðr; f Þ ¼ wT ðr; f Þ R̃con ð f Þwðr; f Þ

ð18Þ

P̂act ðr; n; f Þ ¼ wT ðr; f Þ R̃act ðn; f Þwðr; f Þ

ð19ÞP̂N ðr; n; f Þ ¼ r2 ðn; f ÞwT ðr; f Þwðr; f Þ

ð20Þ

FdB ðr; n; f Þ ¼ 10 log10P̂act ðr; n; f Þ  P̂N ðr; n; f Þ:P̂con ðr; n; f Þ  P̂N ðr; n; f Þð21ÞThis formulation accounts for amplitude differences between
different frequency bands, but its performance may be degraded in
the presence of activity that is more transient. Sources that are
active only briefly may not be adequately captured. Similarly, thespatial filters may not be optimized for sources that change position and orientation over time. Lastly, when analyzing longepochs, this method might be prone to sources active at different
latencies interfering with each other. For example, if one source has
an early response, and another nearby source becomes active laterin the same frequency band, then generating weights from thecovariance of the whole interval may result in degraded reconstruction and poor separation of the two sources.

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Proposed time–frequency optimized beamforming
To overcome the above mentioned limitations, we propose thata custom set of weights w(r, n, f ) be generated from the sample
covariances R̃ act(n, f ) corresponding to each time–frequency window. As in the approaches described above, the data is first passed
through a filter bank and subsequently segmented into overlapping
active windows, τact[n], and control windows, τcon[n]. For optimum time–frequency resolution and beamformer performance, it isdesirable to choose larger time windows for lower frequencies andnarrower time windows for higher frequencies.
R1 ðn; f Þl ðrÞwðr; n; f Þ ¼ T
l ðrÞR1 ðn; f Þl ðrÞ

ð22Þ

where R(n, f ) = [R̃ act(n, f ) + R̃ con( f )] / 2. Then,P̂con ðr; n; f Þ ¼ wT ðr; n; f Þ R̃con ð f Þwðr; n; f Þ

ð23Þ

P̂act ðr; n; f Þ ¼ wT ðr; n; f Þ R̃act ðn; f Þwðr; n; f Þ

ð24ÞP̂N ðr; n; f Þ ¼ r2 ðn; f ÞwT ðr; n; f Þwðr; n; f Þ

ð25Þ

FdB ðr; n; f Þ ¼ 10 log10P̂act ðr; n; f Þ  P̂N ðr; n; f Þ:P̂con ðr; n; f Þ  P̂N ðr; n; f Þð26Þ

Finally, the estimated power of overlapping segments is averaged
to improve numerical stability and better capture transitions insource activity. The procedure is summarized in Fig. 1.
The computational load of the algorithm scales linearly with thenumber of time–frequency bins. In practice, hundreds of weightvectors must be computed to assemble a complete source spectrogram and would require dozens of CPU hours to complete.
However, since the result for each time–frequency window isessentially an independent computation, the time–frequency arrayis well-suited for running on a parallel computing cluster. We used
the shared computing cluster at the California Institute for Quantitative Biomedical Research to generate the results shown in thispaper; the total running time for generating images for all time–
frequency windows is less than 20 min when the cluster is unloaded and all windows can be processed on approximately 300
nodes simultaneously. Upon conclusion of the cluster run, theresults were assembled and visualized with a development versionof our NUTMEG neuromagnetic source reconstruction toolbox(Dalal et al., 2004), freely available from http://bil.ucsf.edu.
The filter bank approach provides an inherent potential advantage over FFT and wavelet-based techniques, since frequency
bins can be of variable size and customized according to theexperimenter's hypothesis. For example, it has been suggested thatthe spectral peak of high gamma activity may vary across subjects
and even within an individual (Crone et al., 2001a; Edwards et al.,2005); therefore, those bands may be defined with a larger
bandwidth. We chose to follow traditional MEG/EEG power band
definitions as best as possible for the experiments presented here:
4–12 Hz (theta–alpha), 12–30 Hz (beta), 30–55 Hz (low gamma),
and 65–90 Hz (high gamma). Additionally, we defined ultrahighfrequency bands at 90–115 Hz, 125–150 Hz, 150–175 Hz, and
185–300 Hz. The power line frequency (60 Hz) and harmonics(120 Hz and 180 Hz) were avoided to reduce noise.
The length of the time windows and width of the frequency bands
must be chosen to ensure that stable and well-conditioned estimates

1689

of the covariance matrices are produced. The parameters used in thispaper were determined with this in mind, but may need to be
adjusted for other MEG systems or data characteristics. See theSupplementary Methods online for an exploration of the effect oftime window length and bandwidth on performance of the proposed
algorithm.
Across-subjects statistics
The significance of activations across subjects was tested withstatistical non-parametric mapping (SnPM) (http://www.sph.umich.
edu/ni-stat/SnPM/). SnPM does not depend on a normal distributionof power change values across subjects and allows correction forfamilywise error of testing at multiple voxels and time–frequency
points. The detailed rationale and procedures of SnPM statistics ofbeamformer images are described elsewhere (Singh et al., 2003). In
short, time–frequency beamformer images for each subject werefirst spatially normalized to the MNI template brain using SPM(http://www.fil.ion.ucl.ac.uk/spm). The three-dimensional averageand variance maps across subjects were calculated at each time–
frequency point. The variance estimates can be noisy for a relatively
low number of subjects, so the variance maps were smoothed with a
20 × 20 × 20 mm3 Gaussian kernel. From this, a pseudo-t statistic was
obtained at each voxel, time window, and frequency band. In
addition, a distribution of pseudo-t statistics was also calculated
from 2N permutations of the original N datasets (subjects). Eachpermutation consisted of two steps: (1) inverting the polarity of thepower change values for some subjects (with 2N possible combinations of negations) and (2) finding the current maximum
pseudo-t value among all voxels and time windows for eachfrequency band. Instead of estimating the significance of each nonpermuted pseudo-t value from an assumed normal distribution, it isthen calculated from the position within the distribution of thesemaximum permuted pseudo-t values. The comparison against
maximum values effectively corrects for the familywise error of
testing multiple voxels and time windows.
Numerical experiments
Data generation
Numerical experiments were conducted to evaluate the proposed
method and compare it with existing methods. The sensor configuration of the 275-channel CTF Omega 2000 biomagnetic measurement system (VSM MedTech, Coquitlam, British Columbia,
Canada) was used. Data were simulated and processed using adevelopment version of NUTMEG (Dalal et al., 2004).
Fifty trials of simulated data were generated, spanning − 750 ms
to 1000 ms per trial, sampled at 1200 Hz. Two 77 Hz sine wave
sources were synthesized and placed at (10, 50, 60) mm and (15,60, 75) mm;1 the phases of each source were assigned randomly
and varied between each other and each trial. The sine waves were
windowed such that they represented ERS activity and were not

1
The numerical experiments used a coordinate system based on a realsubject's head geometry, described as follows: The midpoint between theleft and right preauricular points was defined as the coordinate origin. The
x-axis was directed from the origin through the nasion, while the y-axis was
directed through the left preauricular point and rotated slightly to maintainorthogonality with the x-axis. The z-axis is directed upward perpendicularly
from the xy-plane towards the vertex.

1690

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Fig. 1. Algorithm for optimal time–frequency beamforming. Processing of the combined θ–α band is shown in detail; each of the other frequency bands has asimilar workflow. Note that the algorithm is highly parallel and well-suited to run on high performance computing clusters.

simultaneously active; one source was active from 50 ms to
300 ms, while the other was active from 350 ms to 550 ms. A third
19 Hz source was placed at (25, 30, 100) mm, active from −750 ms
to 50 ms and from 600 ms to 1000 ms to simulate ERD activity.
A sensor lead field was calculated with 5 mm grid spacing
using a single-layer multiple sphere volume conductor as theforward model (Huang et al., 1999) and the Omega 2000's sensor
geometry with respect to a real subject's head-shape. Spontaneous
MEG recordings from a human subject (“brain noise”) were added
to the generated data such that the signal-to-noise (SNR) was equal
to 1. The SNR was defined as the ratio of the Frobenius norm of
the simulated data matrix to that of the brain noise matrix.
Data processing
Covariances for use with the beamformers were generated bycreating a lattice of time–frequency windows. The original data werefirst passed through a bank of 200th–order finite impulse response(FIR) bandpass filters and subsequently split into 29 overlapping
temporal windows with a step size of 25 ms for all bands. In our filter
design, we chose to follow traditional MEG/EEG power band
definitions as best as possible. Theta–alpha band was defined as 4–
12 Hz with 300 ms windows, beta band 12–30 Hz with 200 ms
windows, low gamma 30–55 Hz with 150 ms windows. Additionally, five high gamma bands were defined, avoiding the 60 Hz powerline frequency and its harmonics: 65–90 Hz, 90–115 Hz, 125–150 Hz, 150–175 Hz, 185–300 Hz, all with 100 ms windows.

Finally, sample covariances were calculated for this matrix of time–
frequency windows and averaged over trials:
R̃act ðn; f Þ ¼ hBf ðt act ½nÞBTf ðt act ½nÞi

ð27Þ

Spatial filter weights were computed for each time–frequency
window, and an FdB(r, n, f ) space–time–frequency power map wasassembled as described earlier.
For comparison, the data was processed in three additional ways.
In the first way, which we will term the “broadband” approach, thesimulated data were processed with a conventional minimum
variance beamformer; i.e., a single weight was computed fromunfiltered data using one large active window and a correspondinglarge control window (Eq. (2)). In this case, 0 ms to 500 ms waschosen as the active window, and − 600 ms to − 100 ms was chosen
as the control window. This weight was then applied to the sample
covariances for each time–frequency window to calculate estimatedpower and contrasted with the estimated power of the control period
to generate the final space–time–frequency representation.
The second way was a frequency-dependent beamformer
approach (Eq. (14)). The original simulated data was passed throughthe same filter bank as with our proposed method. However, insteadof segmenting into several time windows, the single large active
window with a corresponding control window was chosen as for thebroadband approach. Thus, weights were computed for filtered data
corresponding to each frequency band. Again, 0 ms to 500 ms was

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

chosen as the active window, with −600 ms to − 100 ms as thecontrol window. Weights, powers, and the final power map weregenerated as with the other two techniques.
Finally, the data was analyzed with sLORETA (Pascual-Marqui,2002) as a representative of minimum norm source reconstruction
techniques. As sLORETA is a nonadaptive spatial filter dependent
only on sensor configuration, the same set of weights was applied tothe sample covariances for each time–frequency window. The
estimated power and contrast with a control period was performed as
described above with the beamformer techniques (Eq. (13)).
Finger movement data
Subjects
Data was collected from 12 right-handed volunteers (6 females
and 6 males, mean age 29.2 years, age range 22–38 years). The
participants were screened for potentially confounding healthconditions and medications. The study protocol was approved by
the UCSF Committee on Human Research.

1691

frontotemporal region (Fig. 2(a)). The electrodes had a 2.3 mm
contact diameter and a center-to-center spacing of 10 mm. Electrodes with an impedance greater than 5 kΩ or exhibiting epileptiform activity were rejected from further analyses. An electrodein the corner of the electrode grid was selected as the reference.
Data was collected with an EEG amplifier (SA Instrumentation,
San Diego, CA) sampling at 2003 Hz with 16-bit resolution. As
with the MEG experiment, patients were asked to move their rightindex finger (RD2) at a self-paced interval of approximately four
seconds for a total of 100 trials. Both patients had corresponding
MEG recordings acquired one day prior to their grid implants. The
recordings were conducted identically as with the healthy volunteers (see above).
Electrodes were localized on individual subject MRIs using
visual identification of landmarks on intraoperative photographs
and backprojection from postimplant X-rays as described by Dalalet al. (submitted for publication) (Fig. 2). Time–frequency analyses

Data acquisition and processing
Data was acquired with a 275-channel CTF Omega 2000 wholehead MEG system from VSM MedTech with a 1200 Hz samplingrate. All post-processing and analysis were performed using adevelopment version of NUTMEG (Dalal et al., 2004). A digital
filter was used to high-pass the data at 1 Hz. Trials containing
eyeblink and movement artifacts were manually rejected.
Subjects were instructed to press the response button with theirright index finger (RD2) at a self-paced interval of approximately
four seconds, acquiring 100 trials. In a subsequent block, thesubjects completed the same task with their left index finger (LD2)instead.
The data was processed as in the above simulation, but with
50 ms window step size due to the length of the epochs. For thebroadband and frequency–domain methods, the active window waschosen to be − 250 ms to 250 ms relative to the button press, with− 950 ms to − 450 ms as the baseline. These windows were chosen
based on typical results in the literature (Pfurtscheller and Neuper,
1992; Jurkiewicz et al., 2006) and our laboratory's extensive
unpublished clinical data.
As with the simulations, a multiple sphere head model was
calculated for each subject at 5 mm resolution based on individual
head shape and relative sensor geometry. Spectral power changes
were statistically tested across subjects with the SnPM method
described above, with p b 0.05 as the threshold for significant
activity.
Intracranial recordings
Preoperative MEG data and corresponding intracranial electrocorticograms (ECoG) were obtained from two patients undergoing
surgical treatment for intractable epilepsy. Intracranial electrodes
were implanted in these patients for preresection seizure localization and functional mapping of critical language and motor areas.
The study protocol, approved by the UCSF and UC Berkeley
Committees on Human Research, did not interfere with the ECoG
recordings made for clinical purposes and presented minimal risk
to the subjects. Upon informed consent, the experiments wereconducted while the patient was alert and on minimal medication.
The implants consisted of an 8 × 8 grid of platinum–iridium
electrodes (Ad-Tech Medical, Racine, WI) placed over the left

Fig. 2. (a) Example of a typical frontotemporal ECoG montage in an
intractable epilepsy patient. The implant consists of an 8 × 8 electrode gridwith 10 mm center-to-center spacing between electrodes. (b) Lateral X-ray
radiograph of the same patient showing electrode locations. The surgicalphotograph was used to annotate the locations of visible electrodes on an
MRI rendering, while the coordinates of hidden electrodes were found using
X-ray backprojection to the MRI-derived brain surface (Dalal et al.,submitted for publication).

1692

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Fig. 3. At top is the spectrogram corresponding to the three simulated sources. In the rows below are the reconstruction results using sLORETA, the broadband
beamformer, the frequency domain beamformer, and the proposed time–frequency beamformer. In each of those panels, the crosshairs mark the spatiotemporalpeak for the reconstructed source, with the corresponding spectrogram shown below it. The time–frequency window plotted on the MRI is highlighted on the
spectrogram. The functional maps are thresholded at 50% of the maximum power (in dB) for the beamformer variants and 75% for sLORETA.

Fig. 4. Shown at top are the grand average reconstruction results for right index finger movement using the broadband beamformer, the frequency domain
beamformer, and the proposed time–frequency beamformer. The functional maps are superimposed on the MNI template brain and are statistically thresholded at
p b 0.05 (corrected). In each of the panels, the crosshairs mark the spatiotemporal peak for the reconstructed source, with the corresponding spectrogram shown
below it. The functional map plotted on the MRI corresponds to the time–frequency window highlighted on the spectrogram. Note that the frequency–domain
beamformer localized peaks similar to the other methods, but grossly overestimated the statistically significant spatial extent of the late beta ERS, likely due tothe large baseline shift of inactive voxels.

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

1693

1694

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

of ECoG data were performed using the event-related spectral
perturbation (ERSP) method (Makeig, 1993). Time courses for the
power of single trial data were generated for each frequency bandusing a Gaussian filter bank and the Hilbert transform (Edwards,
2007); after averaging across trials, the power time courses weredivided by the mean baseline spectrum to generate the ERSP. These
results were converted to decibels and then rebinned into the sametime–frequency windows used to analyze the MEG data for ease ofcomparison.
Results
Numerical experiments
The sLORETA method produced relatively blurry results for all
three simulated sources, with peaks on the periphery of the definedvolume of interest in each case (see Fig. 3). The reconstructions
were not of sufficiently high fidelity to appreciably distinguish thespectrograms of the different sources. Several regularization parameters were tested with similar results.
The broadband beamformer correctly placed the peak of beta
ERD at (25, 30, 100) mm (see Fig. 3). However, the spatiotemporalextent of both high gamma ERS sources were not as cleanly
resolved. The first source was placed at (20, 50, 60) mm peaking
over 150–175 ms, while the second source was placed at (20, 55,70) mm, peaking over 450–475 ms. Additionally, the spatial extent
of all sources was blurred.
The frequency domain beamformer found the correct location of
the beta ERD, resolving a more focal peak than the broadband
beamformer (see Fig. 3). It also found the correct locations for both
high gamma ERS sources. However, the activation was spatiallyblurred and attenuated for the high gamma ERS sources, especiallyover 300–350 ms when one source tapers off and other tapers on.
Additionally, the spatiotemporal extent of all three sources wascompromised. The spectrogram computed for (15, 60, 75) mm shows
contamination from the (10, 50, 60) mm source and vice versa.
Finally, we applied our proposed technique to the data (see Fig.
3). As expected, the beta ERD was accurately resolved. Both highgamma ERS sources were accurately localized and their temporalextents accurately captured. Virtually no contaminations between
the two source locations were observed on their respective
spectrograms. This method provided the best reconstruction of thesimulated data.
Finger movement data
The characteristic beta band power decrease in contralateral
sensorimotor cortex was observed and reached statisticalsignificance across subjects for all three beamformer variants(see Fig. 4 for corresponding MNI coordinates and corrected pvalues for right index finger movement). However, only lowamplitude early time windows near − 400 ms were significant
for the broadband beamformer. In contrast, significant contralateral activation was observed over − 500 ms to 250 ms withboth the proposed time–frequency beamformer and thefrequency–domain beamformer, although results were more
spatially focal for the proposed method. Additionally, both ofthese methods revealed significant beta band power decreases in
ipsilateral sensorimotor cortex and ipsilateral secondary somatosensory cortex approximately 0 ms to 200 ms after movement
onset.

The contralateral decrease in beta power was followed by asignificant contralateral beta rebound for all three methods. Again,
the time–frequency beamformer performed the best, with arelatively focal activation area. The broadband beamformer
revealed a peak in sensorimotor cortex, but the spatial extent of
the activation extended into areas both implausibly deep as well as
outside the brain. The frequency–domain beamformer placed thepeak nearby, but grossly overestimated the statistically significant
spatial extent, apparently due to a large baseline shift evident invoxels distant from motor cortex. The time–frequency beamformer
depicts a relatively focal activation in contralateral sensorimotor
cortex. (Individual results for many subjects also showed an
ipsilateral beta rebound, but this did not reach statisticalsignificance across subjects.) It also found an increase in beta
power peaking at (5, − 5, 65) mm (MNI coordinates, p b 0.038,
corrected), corresponding to activation of the supplementary motor
area (SMA) (not shown).
Interestingly, both the frequency–domain beamformer and thetime–frequency beamformer localized a focal, statistically significant high gamma (65–90 Hz) peak in sensorimotor cortex. This
activity was found to be more spatially focal and temporally boundto the movement. No significant high gamma activity was observed with the broadband beamformer.
Similarly, the proposed technique revealed similar activity forleft index finger movement (Fig. 5). The typical beta banddesynchronization and late rebound as well as high gamma activity
were found in right sensorimotor cortex, reaching statisticalsignificance across subjects.
Activation of the cerebellum was also found in 9 of 12 healthy
volunteers and in both of the patients (see Fig. 6). While thespatiotemporal extent and particular frequency content of cerebellar activations exhibited considerable variability across subjects
and did not reach statistical significance in our across-subject
analyses with whole-brain multiple comparison correction, we didobserve that our method found consistent high-frequency sourcesin the cerebellum in either the 65–90 Hz or 90–115 Hz bands.
Examples of distinct cerebellar responses from two subjects areshown in Fig. 6; see Fig. 7 for responses from the two patients.
Intracranial recordings
As shown in Fig. 7, several locations showing ECoG activity
during the right finger movement task were also found with theproposed MEG time–frequency beamformer method and exhibited
fairly similar spectrogram patterns. Table 1 lists the coordinates ofeach peak in the grid coverage area for beta and high gamma
activations for both patients. MEG peaks were found between
2.8 mm and 10.4 mm from eight ECoG peaks, while two adjacentelectrodes showing low-amplitude beta ERD and one electrodeshowing high gamma ERS did not have corresponding MEG
activations.
Note that the MEG reconstruction for both patients show thelargest-amplitude beta desynchronization and high gamma synchronization in left primary motor cortex and the cerebellum inaccordance with the across-subjects analyses above, but these areas
were not covered by the ECoG grid in either patient; therefore, theECoG analyses show only lower-amplitude secondary areas ofactivation which tend to result in blurrier MEG activations.
Nevertheless, the ECoG analyses supported the validity of MEG
reconstructions of these secondary activations, taking into accountthe 1 cm spacing and cortical surface placement of the ECoG grid

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Fig. 5. Shown above are the grand average reconstruction results for left index finger movement using the proposed time–frequency beamformer, superimposed on the MNI template brain. The functional maps aresuperimposed on the MNI template brain and are statistically thresholded at p b 0.05 (corrected). In each panel, the crosshairs mark the spatiotemporal peak for the reconstructed source, with the corresponding
spectrogram shown below it. The functional map plotted on the MRI corresponds to the time–frequency window highlighted on the spectrogram.

1695

1696

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Fig. 6. Above, examples of cerebellum activation for finger movement in two subjects. Above left are the results for RD2 movement in one subject. Above right
are the results for LD2 movement in a different subject. Both functional maps are thresholded at 75% of the maximum power (in dB).

as well as spatiotemporal blurring inherent to the beamformer
technique.
Discussion
We have shown that, with our novel time–frequency optimized
beamformer techniques, MEG can resolve sources of transient
power changes across multiple frequency bands, including highgamma activity. The method was validated with across-subjects
statistics and intracranial recordings.
Some secondary activity revealed by the ECoG analyses was notobserved with the MEG source reconstructions; these sources may
have activated a small cortical region and/or were not optimallyoriented for detection by MEG sensor arrays. Additionally, MEG
source reconstructions for any given voxel are linear combinationsof activity from multiple nearby sources due to spatiotemporal blurand may explain minor spectrogram differences as compared to
ECoG. The degree of spatial blur depends on various factors,
especially SNR as well as the true spatial extent of the sources.
Adaptive spatial filter weights computed in the traditionalmanner from unfiltered or wideband data are inherently biased
towards resolving low-frequency brain activity due to the power lawof typical electrophysiological data. By creating a set of weights
customized for each time–frequency window, higher frequency
sources may be characterized with much greater fidelity. Additionally, segmenting the data into time windows can better capture thetemporal extent of oscillatory modulations as well as allow forsources to change position and orientation. This is particularly
important for experiment designs with long interstimulus intervals
that yield several hundred milliseconds of data per epoch.
In using our proposed method, care must be taken to choose
window lengths and bandwidths that ensure stable and wellconditioned covariance estimates. Optimal parameters will varyconsiderably depending on the type of experiment, backgroundnoise level, number of sensors, and number of trials. Too few timesamples or too little bandwidth would result in poor covarianceestimates and severely degrade performance of the beamformer
reconstruction (Brookes et al., 2008). As we demonstrated in the
Supplementary Methods section (available online), a compromisemust be made between bandwidth and time window length. The
ultimate parameter choice, then, must be driven by experimental
hypotheses. It must be considered that real sources are unlikely tostay active for several hundreds of milliseconds at a time, making

extremely narrow bandwidths impractical. Conversely, a source at
a given location may generate a power increase in one band
simultaneously with a power decrease in another band (as with
common gamma power increase commonly observed in tandemwith beta power decrease); performance would suffer if both eventswere contained by a single wide frequency band. Finally, total data
length (the product of time window length with number of trials)
had a direct impact on performance. Therefore, experiments should
be designed with the duration of hypothesized activations andbands of interest in mind, increasing the number of trials acquired
as necessary.
Noise is known to significantly impact the performance ofminimum-norm-based methods by increasing localization bias anddecreasing spatial resolution (Greenblatt et al., 2005; Sekihara etal., 2005), and this likely explains the presented sLORETA
results; ultimately, the regularization parameters and method are
critical in the presence of noise. Perhaps a similar approach to the
proposed method can be taken; i.e., regularization parameters canbe customized for different time–frequency segments, creating a
hybrid adaptive–nonadaptive source reconstruction technique.
This requires additional investigations beyond the scope of thispaper.
ECoG has been shown to clearly resolve high gamma (N 60 Hz)
activity and suggests it is more spatiotemporally focal than lowerfrequency activity (Crone et al., 1998, 2001a,b; Edwards et al.,2005; Canolty et al., 2006). Recently, high gamma activity has
been gaining attention in the MEG/EEG literature as well (Kaiser
et al., 2002; Hoogenboom et al., 2006; Vidal et al., 2006; Siegel
et al., 2007; Osipova et al., 2006). While increases in high gamma
power may coincide with decreases in beta power in many cases,
high gamma may be a better indicator of task-specific neuralprocessing in local cortical circuits since it is found to be more
focused spatially and temporally. The hand motor data we present
in this paper supports this hypothesis. Additionally, many studies
have recently shown that high gamma activity is positively
correlated with the hemodynamic response measured by functional
MRI (fMRI) (Logothetis et al., 2001; Mukamel et al., 2005;
Niessing et al., 2005; Brovelli et al., 2005; Hoogenboom et al.,2006; Lachaux et al., 2007). Finally, higher-frequency bands may
be less likely to be temporally correlated even if they are
simultaneously active, and may thereby naturally circumvent the
known limitation of beamformer techniques to resolve highly
temporally correlated sources (Sekihara et al., 2002).

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700
Fig. 7. Shown above are the right finger (RD2) movement activity for two intractable epilepsy patients, using both time–frequency analyses from an 8 × 8 intracranial electrode grid and the corresponding results frompreoperative magnetoencephalography with the proposed time–frequency beamformer. The spectrogram corresponds to the circled spatial location, while the functional maps show the spatial extent of activation for
the indicated time window and frequency band. The orange outline indicates the region covered by the intracranial electrode grid. Note that MEG reveals strong primary motor cortex and cerebellum activity, but theseareas were not covered with electrodes in either patient; instead, lower-amplitude secondary activations are compared between the two methods.
1697

1698

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Table 1
ECoG peaks vs. MEG reconstruction peaks
Patient Band

ECoG coordinates MEG coordinates
(mm)
(mm)

1
1
1
1
1
1
1
2
2
2
2

−50.6, 3.8, 28.2
−55.2, 13.2, 9.4−54.9, − 0.7, 20.6
−57.7, − 5.1, 0.0−57.7, 0.6, 7.7−54.9, − 0.7, 20.6
−53.7, 20.3, 3.9−56.9, − 24.4, 30.6
−58.3, 2.0, 14.1
−55.0, − 33.7, 33.3
−64.4, − 26.4, 23.5

beta
beta
beta
beta
beta
high gamma
high gamma
beta
beta
beta
high gamma

Difference(mm)− 45.8, −1.2, 25.5
7.5
− 50.3, 13.8, 15.3
7.7
− 50.8, −1.1, 25.4
6.3
–
–
–
–
− 51.0, −3.1, 30.0 10.4––
− 51.3, − 28.9, 36.7 9.5− 61.0, 0.3, 8.9
6.1− 51.2, − 27.5, 32.0 7.4− 66.2, − 24.5, 22.5 2.8

ECoG electrode locations with activity in the beta and high gamma bands arelisted along with the nearest peaks found from the MEG time-frequency
beamformer reconstruction. Note that coordinates given are in each patient's
native MRI space (rather than MNI coordinates) in order to accurately
characterize Euclidean distances between ECoG and MEG peaks, given in
the last column.

Other ECoG studies also show motor ERS in bands greater than
60 Hz (Ohara et al., 2000; Pfurtscheller et al., 2003) and even up to
200 Hz (Leuthardt et al., 2004; Brovelli et al., 2005; Crone et al.,2006) in the same region we observed with our MEG technique.
Additionally, the postmovement beta rebound has been observed inboth ECoG (Pfurtscheller et al., 1996; Sochůrková et al., 2006) andMEG (Jurkiewicz et al., 2006).
Our method suggested activations in the cerebellum for most of
the healthy subjects and both epilepsy patients, though it did notreach statistical significance across subjects, likely due to individual
variability in precise location, latency, and frequency. We speculatethat both the sensor configuration and existing head models are not
optimized for accuracy in the cerebellar region. Currently available
MEG sensor arrays may not provide adequate coverage that far
down the head with normal subject positioning. Furthermore,evidence from fMRI studies (Grodd et al., 2001; Hülsmann et al.,
2003; Thickbroom et al., 2003; Dhamala et al., 2003; Dimitrova
et al., 2006) suggests that the anterior cerebellum may be the mostactive, placing the neural generators fairly distant from the sensors
and significantly lowering the SNR of the signals. Additionally, thestrategy employed by individual subjects in pacing their finger
movements may have introduced variability in the quality and extent
of activation due to the cerebellum's role in timing and rhythm (Ivryand Keele, 1989; Dhamala et al., 2003; Lotze et al., 2003). Finally,
existing MEG/EEG head models focus on cerebral hemispheres anddo not explicitly account for the structure of the cerebellum or its rolein generating signals. As such, they may introduce large lead field
inaccuracies in the region of the cerebellum, severely degrading theperformance of spatial filtering techniques. Perhaps more sophisticated models based on boundary element modeling (BEM) or finite
element modeling (FEM) are needed to improve fidelity in the
cerebellum and other deep brain structures.
Previous MEG/EEG studies have suggested coherence between
the cerebellum and cerebral cortex in the alpha and beta bands (Grosset al., 2002; Pollok et al., 2005). However, the activations found in
this study suggest that the cerebellum may exhibit oscillatory activity
at much higher frequencies that are not necessarily coherent withother locations, in accordance with speculation by Niedermeyer(2004) and the classic experiments of Adrian (1935), Dow (1938),

Ten Cate and Wiggers (1942), and Pellet (1967). The use of space–time–frequency methods for analyzing MEG/EEG data may finallyallow the cerebellum's electrical activity to be independently studied
noninvasively.
In addition to the hand motor data presented in this paper, other
MEG studies by our group show that our method can reveal more
complex cognitive processes related to learning, decision-making,and memory (Dalal et al., 2005; van Wassenhove and Nagarajan,
2006, 2007; Hinkley, 2007; Guggisberg et al., 2008).
The technique we propose can be customized according to thepreferences of the experimenter. For example, the frequency bands
and time windows can be adjusted depending on the expected SNR
and trial-to-trial variability of the experiment. Any typical filter
type can be used to construct the filter banks; an experimenter may
prefer to substitute filters with different properties than we havechosen or even wavelet-based filters. Finally, since the power of
the active windows, control window, and noise are preserved in the
final results, the contrast type may be selected by the end user.
Rather than an F-ratio contrast, a t-test (difference) or theuncontrasted power time course may be selected instead.
This type of analysis does yield a large amount of information—atime–frequency spectrogram for every spatial location implies fivedimensions of output data! Therefore, we have implemented aninteractive time–frequency viewer into our software package
NUTMEG to help make navigation of the results more intuitive.
Future directions may include developing factor analysis techniques
to help mine the rich output afforded by five-dimensional space–time–frequency analyses.
Acknowledgments
The authors would like to thank Johanna M. Zumer, Susanne
M. Honma, Virginie van Wassenhove, Leighton B.N. Hinkley,
John F. Houde, Julia P. Owen, Maryam Soltani, and Mary M.
Mantle for invaluable assistance and feedback, as well as Jeff
Block and Jason Crane for critical advice on the use of parallel
computing resources. SSD was supported in part by NIH grant F31
DC006762 and SSN was supported in part by NIH grants R01
DC004855 and DC006435.
Appendix A. Supplementary data
Supplementary data associated with this article can be found, inthe online version, at doi:10.1016/j.neuroimage.2008.01.023.

References
Adrian, E.D., 1935. Discharge frequencies in the cerebral and cerebellar
cortex. J. Physiol. 83, 32P–33P.
Berger, H., 1930. Über das elektrenkephalogramm des menschen. J. Psychol.
Neurol. 40, 160–179.
Brookes, M.J., Gibson, A.M., Hall, S.D., Furlong, P.L., Barnes, G.R.,
Hillebrand, A., Singh, K.D., Holliday, I.E., Francis, S.T., Morris, P.G.,
2004. A general linear model for MEG beamformer imaging. NeuroImage 23, 936–946.
Brookes, M.J., Vrba, J., Robinson, S.E., Stevenson, C.M., Peters, A.M.,
Barnes, G.R., Hillebrand, A., Morris, P.G., 2008. Optimising experimental
design for MEG beamformer imaging. NeuroImage 39, 1788–1802.
Brovelli, A., Lachaux, J.-P., Kahane, P., Boussaoud, D., 2005. High gamma
frequency oscillatory activity dissociates attention from intention in thehuman premotor cortex. NeuroImage 28, 154–164.

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700
Canolty, R.T., Edwards, E., Dalal, S.S., Soltani, M., Nagarajan, S.S., Kirsch,
H.E., Berger, M.S., Barbaro, N.M., Knight, R.T., 2006. High gamma
power is phase-locked to theta oscillations in human neocortex. Science
313, 1626–1628.
Cheyne, D., Gaetz, W., Garnero, L., Lachaux, J.-P., Ducorps, A., Schwartz, D.,
Varela, F.J., 2003. Neuromagnetic imaging of cortical oscillations accompanying tactile stimulation. Brain Res. Cogn. Brain Res. 17, 599–611.
Crone, N.E., Miglioretti, D.L., Gordon, B., Lesser, R.P., 1998. Functional
mapping of human sensorimotor cortex with electrocorticographic
spectral analysis: II. Event-related synchronization in the gamma band.
Brain 121, 2301–2315.
Crone, N.E., Boatman, D., Gordon, B., Hao, L., 2001a. Induced
electrocorticographic gamma activity during auditory perception. Clin.
Neurophysiol. 112, 565–582.
Crone, N.E., Hao, L., Hart, J., Boatman, D., Lesser, R.P., Irizarry, R.,
Gordon, B., 2001b. Electrocorticographic gamma activity during word
production in spoken and sign language. Neurology 57, 2045–2053.
Crone, N.E., Sinai, A., Korzeniewska, A., 2006. Chapter 19: Highfrequency gamma oscillations and human brain mapping with electrocorticography. Prog. Brain Res. 159, 275–295.
Dalal, S.S., Zumer, J.M., Agrawal, V., Hild, K.E., Sekihara, K., Nagarajan,
S.S., 2004. NUTMEG: a neuromagnetic source reconstruction toolbox.
Neurol. Clin. Neurophysiol. 52.
Dalal, S.S., Edwards, E., Kirsch, H.E., Canolty, R.T., Soltani, M., Barbaro,
N.M., Knight, R.T., Nagarajan, S.S., 2005. Spatiotemporal dynamics ofcortical networks preceding finger movement and speech production[abstract]. Proceedings of the 16th Meeting of the International Society
for Brain Electromagnetic Topography. Brain Topogr. 18, 128.
Dalal, S.S., Edwards, E., Kirsch, H.E., Barbaro, N.M., Knight, R.T.,
Nagarajan, S.S., submitted for publication. Localization of neurosurgicallyimplanted electrodes via photograph–MRI-radiograph coregistration.
Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R.M., 2004.
Mapping human brain function with MEG and EEG: methods andvalidation. NeuroImage 23 (Suppl 1), S289–S299.
Dhamala, M., Pagnoni, G., Wiesenfeld, K., Zink, C.F., Martin, M., Berns, G.S.,
2003. Neural correlates of the complexity of rhythmic finger tapping.
NeuroImage 20, 918–926.
Dimitrova, A., de Greiff, A., Schoch, B., Gerwig, M., Frings, M., Gizewski,
E.R., Timmann, D., 2006. Activation of cerebellar nuclei comparing
finger, foot and tongue movements as revealed by fMRI. Brain Res.
Bull. 71, 233–241.
Dow, R.S., 1938. The electrical activity of the cerebellum and its functionalsignificance. J. Physiol. 94, 67–86.
Edwards, E., May 2007. Electrocortical activation and human brain
mapping. Ph.D. thesis, University of California, Berkeley.
Edwards, E., Soltani, M., Deouell, L.Y., Berger, M.S., Knight, R.T., 2005.
High gamma activity in response to deviant auditory stimuli recorded
directly from human cortex. J. Neurophysiol. 94, 4269–4280.
Gaetz, W., Cheyne, D., 2006. Localization of sensorimotor cortical rhythms
induced by tactile stimulation using spatially filtered MEG. NeuroImage
30, 899–908.
Greenblatt, R.E., Ossadtchi, A., Pflieger, M.E., 2005. Local linear estimators
for the bioelectromagnetic inverse problem. IEEE Trans. Signal Process.
53, 3403–3412.
Grodd, W., Hülsmann, E., Lotze, M., Wildgruber, D., Erb, M., 2001.
Sensorimotor mapping of the human cerebellum: fMRI evidence ofsomatotopic organization. Hum. Brain Mapp. 13, 55–73.
Gross, J., Kujala, J., Hamalainen, M., Timmermann, L., Schnitzler, A.,
Salmelin, R., 2001. Dynamic imaging of coherent sources: studying
neural interactions in the human brain. Proc. Natl. Acad. Sci. U. S. A. 98,
694–699.
Gross, J., Timmermann, L., Kujala, J., Dirks, M., Schmitz, F., Salmelin, R.,
Schnitzler, A., 2002. The neural basis of intermittent motor control inhumans. Proc. Natl. Acad. Sci. U. S. A. 99, 2299–2302.
Guggisberg, A.G., Dalal, S.S., Findlay, A.M., Nagarajan, S.S., 2008. Highfrequency oscillations in distributed neural networks reveal the
dynamics of human decision making. Front. Hum. Neurosci. 1.

1699

Hämäläinen, M., Hari, R., IImoniemi, R.J., Knuutila, J., Lounasmaa, O.V.,
1993. Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod.
Phys. 65, 413–497.
Hanslmayr, S., Klimesch, W., Sauseng, P., Gruber, W., Doppelmayr, M.,
Freunberger, R., Pecherstorfer, T., Birbaumer, N., 2007. Alpha phase
reset contributes to the generation of ERPs. Cereb. Cortex 17 (1), 1–8.
Hinkley, L.B.N., September 2007. The role of human posterior parietal
cortex in sensory, motor and cognitive function. Ph.D. thesis, Universityof California, Davis.
Hoogenboom, N., Schoffelen, J.-M., Oostenveld, R., Parkes, L.M., Fries, P.,
2006. Localizing human visual gamma-band activity in frequency, timeand space. NeuroImage 29, 764–773.
Huang, M.X., Mosher, J.C., Leahy, R.M., 1999. A sensor-weighted
overlapping-sphere head model and exhaustive head model comparison
for MEG. Phys. Med. Biol. 44, 423–440.
Hülsmann, E., Erb, M., Grodd, W., 2003. From will to action: sequential cerebellar contributions to voluntary movement. NeuroImage 20, 1485–1492.
Ivry, R.B., Keele, S.W., 1989. Timing functions of the cerebellum. J. Cogn.
Neurosci. 1, 136–152.
Jensen, O., Vanni, S., 2002. A new method to identify multiple sources ofoscillatory activity from magnetoencephalographic data. NeuroImage
15, 568–574.
Jurkiewicz, M.T., Gaetz, W.C., Bostan, A.C., Cheyne, D., 2006. Postmovement beta rebound is generated in motor cortex: evidence fromneuromagnetic recordings. NeuroImage 32, 1281–1289.
Kaiser, J., Lutzenberger, W., Ackermann, H., Birbaumer, N., 2002.
Dynamics of gamma-band activity induced by auditory pattern changes
in humans. Cereb. Cortex 12, 212–221.
Lachaux, J.-P., Fonlupt, P., Kahane, P., Minotti, L., Hoffmann, D., Bertrand,O., Baciu, M., 2007. Relationship between task-related gamma
oscillations and BOLD signal: new insights from combined fMRI and
intracranial EEG. Hum. Brain. Mapp. 28, 1368–1375.
Leuthardt, E.C., Schalk, G., Wolpaw, J.R., Ojemann, J.G., Moran, D.W.,
2004. A brain–computer interface using electrocorticographic signals inhumans. J. Neural Eng. 1, 63–71.
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.
Lotze, M., Scheler, G., Tan, H.-R.M., Braun, C., Birbaumer, N., 2003. The
musician's brain: Functional imaging of amateurs and professionals
during performance and imagery. NeuroImage 20, 1817–1829.
Makeig, S., 1993. Auditory event-related dynamics of the EEG spectrum
and effects of exposure to tones. Electroencephalogr. Clin. Neurophysiol. 86, 283–293.
Michalewski, H.J., Prasher, D.K., Starr, A., 1986. Latency variability andtemporal interrelationships of the auditory event-related potentials (N1,P2, N2, and P3) in normal subjects. Electroencephalogr. Clin. Neurophysiol. 65, 59–71.
Mukamel, R., Gelbard, H., Arieli, A., Hasson, U., Fried, I., Malach, R.,
2005. Coupling between neuronal firing, field potentials, and FMRI in
human auditory cortex. Science 309, 951–954.
Niedermeyer, E., 2004. The electrocerebellogram. Clin. EEG Neurosci. 35,
112–115.
Niessing, J., Ebisch, B., Schmidt, K.E., Niessing, M., Singer, W., Galuske,
R.A.W., 2005. Hemodynamic signals correlate tightly with synchronizedgamma oscillations. Science 309, 948–951.
Ohara, S., Ikeda, A., Kunieda, T., Yazawa, S., Baba, K., Nagamine, T., Taki,
W., Hashimoto, N., Mihara, T., Shibasaki, H., 2000. Movement-relatedchange of electrocorticographic activity in human supplementary motor
area proper. Brain 123, 1203–1215.
Osipova, D., Takashima, A., Oostenveld, R., Fernández, G., Maris, E.,
Jensen, O., 2006. Theta and gamma oscillations predict encoding andretrieval of declarative memory. J. Neurosci. 26, 7523–7531.
Pascual-Marqui, R.D., 2002. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find.
Exp. Clin. Pharmacol. 24 (Suppl D), 5–12.

1700

S.S. Dalal et al. / NeuroImage 40 (2008) 1686–1700

Pellet, J., 1967. L'électrocérébellogramme vermien au cours des états deveille et de sommeil. Brain Res. 5, 266–270.
Pfurtscheller, G., Aranibar, A., 1977. Event-related cortical desynchronization detected by power measurements of scalp EEG. Electroencephalogr.
Clin. Neurophysiol. 42, 817–826.
Pfurtscheller, G., Neuper, C., 1992. Simultaneous EEG 10 Hz desynchronization and 40 Hz synchronization during finger movements.
NeuroReport 3, 1057–1060.
Pfurtscheller, G., Stancák, A., Neuper, C., 1996. Post-movement beta
synchronization. A correlate of an idling motor area? Electroencephalogr. Clin. Neurophysiol. 98, 281–293.
Pfurtscheller, G., Graimann, B., Huggins, J.E., Levine, S.P., Schuh, L.A.,
2003. Spatiotemporal patterns of beta desynchronization and gammasynchronization in corticographic data during self-paced movement.
Clin. Neurophysiol. 114, 1226–1236.
Pollok, B., Gross, J., Müller, K., Aschersleben, G., Schnitzler, A., 2005. The
cerebral oscillatory network associated with auditorily paced finger
movements. NeuroImage 24, 646–655.
Robinson, S.E., Vrba, J., 1999. Functional neuroimaging by synthetic
aperture magnetometry. In: Yoshimoto, T., Kotani, M., Kuriki, S.,
Karibe, H., Nakasato, N. (Eds.), Recent Advances in Biomagnetism.
Tohoku University Press, Sendai, pp. 302–305.
Sekihara, K., Scholz, B., 1996. Generalized Wiener estimation of threedimensional current distribution from biomagnetic measurements. In:
Aine, C.J. (Ed.), Biomag 96: Proceedings of the Tenth International
Conference on Biomagnetism. Springer-Verlag, pp. 338–341.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., Miyashita, Y., 2001.
Reconstructing spatio-temporal activities of neural sources using an MEG
vector beamformer technique. IEEE Trans. Biomed. Eng. 48, 760–771.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., 2002. Performanceof an MEG adaptive-beamformer technique in the presence of correlated
neural activities: effects on signal intensity and time–course estimates.
IEEE Trans. Biomed. Eng. 49, 1534–1546.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., 2004. Asymptotic
SNR of scalar and vector minimum-variance beamformers forneuromagnetic source reconstruction. IEEE Trans. Biomed. Eng. 51,
1726–1734.
Sekihara, K., Sahani, M., Nagarajan, S.S., 2005. Localization bias andspatial resolution of adaptive and non-adaptive spatial filters for MEG
source reconstruction. NeuroImage 25, 1056–1067.
Siegel, M., Donner, T.H., Oostenveld, R., Fries, P., Engel, A.K., 2007. Highfrequency activity in human visual cortex is modulated by visual motion
strength. Cereb Cortex 17, 732–741.
Singh, K.D., Barnes, G.R., Hillebrand, A., Forde, E.M.E., Williams, A.L.,
2002. Task-related changes in cortical synchronization are spatiallycoincident with the hemodynamic response. NeuroImage 16,
103–114.
Singh, K.D., Barnes, G.R., Hillebrand, A., 2003. Group imaging of taskrelated changes in cortical synchronisation using nonparametricpermutation testing. NeuroImage 19, 1589–1601.
Sochůrková, D., Rektor, I., Jurák, P., Stančák, A., 2006. Intracerebral
recording of cortical activity related to self-paced voluntary movements:
a Bere itschaftspotential and event-related desynchronization/synchronization. SEEG study. Exp. Brain Res. 173, 637–649.
Ten Cate, J., Wiggers, N., 1942. On the occurrence of slow waves in theelectrocerebellogram. Arch. Neerl. Physiol. l'Homme Anim. 26, 433–435.
Thickbroom, G.W., Byrnes, M.L., Mastaglia, F.L., 2003. Dual representation of the hand in the cerebellum: activation with voluntary and passive
finger movement. NeuroImage 18, 670–674.
Van Veen, B.D., Buckley, K.M., April 1988. Beamforming: a versatile
approach to spatial filtering. IEEE ASSP Mag. 5, 4–24.
van Wassenhove, V., Nagarajan, S.S., April 2006. Auditory–visual habituationin perceptual recalibration [abstract]. Cognitive Neuroscience Society
Annual Meeting Program. J. Cogn. Neurosci. 18 (Suppl 1), 594.
van Wassenhove, V., Nagarajan, S.S., 2007. Auditory cortical plasticity inlearning to discriminate modulation rate. J. Neurosci. 27, 2663–2672.
Vidal, J.R., Chaumon, M., O'Regan, J.K., Tallon-Baudry, C., 2006. Visualgrouping and the focusing of attention induce gamma-band oscillations at
different frequencies in human magnetoencephalogram signals. J. Cogn.
Neurosci. 18, 1850–1862.