www.elsevier.com/locate/ynimg
NeuroImage 30 (2006) 400 – 416

A graphical model for estimating stimulus-evoked brain responses
from magnetoencephalography data with large backgroundbrain activity
Srikantan S. Nagarajan,a,* Hagai T. Attias,b Kenneth E. Hild II,a and Kensuke Sekihara caBiomagnetic Imaging Laboratory, Department of Radiology, University of California at San Francisco, San Francisco, CA 94122, USA
Golden Metallic Inc., San Francisco, CA 94147, USAc
Department of Systems Design and Engineering, Tokyo Metropolitan University, Asahigaoka 6-6, Hino, Tokyo 191-0065, Japanb

Received 16 April 2005; revised 27 September 2005; accepted 29 September 2005
Available online 19 December 2005

This paper formulates a novel probabilistic graphical model for noisy
stimulus-evoked MEG and EEG sensor data obtained in the presence oflarge background brain activity. The model describes the observed datain terms of unobserved evoked and background factors with additive
sensor noise. We present an expectation maximization (EM) algorithm
that estimates the model parameters from data. Using the model, thealgorithm cleans the stimulus-evoked data by removing interference
from background factors and noise artifacts and separates those data
into contributions from independent factors. We demonstrate on realand simulated data that the algorithm outperforms benchmark methods
for denoising and separation. We also show that the algorithm improves
the performance of localization with beamforming algorithms.
D 2005 Elsevier Inc. All rights reserved.
Keywords: Magnetoencephalography; Electroencephalography; Graphicalmodels; EM algorithm

Introduction
Electromagnetic source imaging (ESI), the reconstruction of thespatiotemporal activation of brain sources from MEG and EEG
data, is increasingly being used in numerous studies of humancognition in both normals and in various clinical populations(Hmlinen et al., 1993; Baillet et al., 2001; Ebersole, 1997). The
major advantage of ESI over other noninvasive functional brain
imaging techniques is the ability to obtain valuable information
about neural dynamics with high temporal resolution on the orderof milliseconds. A major problem in ESI is that MEG and EEG
measurements, which use sensors located outside the brain,
generally contain not only signals associated with brain sourcesof interest, but also signals from other sources such as spontaneous* Corresponding author.
E-mail address: sri@radiology.ucsf.edu (S.S. Nagarajan).
Available online on ScienceDirect (www.sciencedirect.com).
1053-8119/$ - see front matter D 2005 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2005.09.055

brain activity, eye blinks and other biological and non-biologicalsources of artifacts. Interference signals from these sources overlap
spatially and temporally with those from the brain sources ofinterest, making it difficult to obtain accurate reconstructions.
Many approaches have been taken to address the problems of
reducing noise from artifacts and from interference signals due tobackground brain activity, with varying degrees of success. First,
event-related averaging of multiple trials may be used to reduce thecontributions from non-stimulus-locked and non-phase-locked
activity. Averaging reduces the contribution of background brain
activity on the order of the square root of the number of trials. This
is often performed after rejection of manually selected trials that
contain conspicuous artifacts, such as eye blinks or saccades.
However, averaging and manual trial segment rejection require thatthe data be collected for a large number of trials, which effectively
limits the number of stimulus conditions in a particular experiment.
Second, model-based approaches have been used to either
model the neural generator (de Munck et al., 2004; Bijma et al.,2003, 2004; Baryshnikov et al., 1981) or to model artifacts (Berg etal., 1991). The former involves defining a parametric model for the
event-related brain response and the background noise andestimating parameters using a maximum-likelihood methods. The
latter involves denoising of the data by identifying and removing
portions of data that fit to previously estimated models of an arrayof artifacts. With the latter technique, only irregularities that havebeen previously identified and modeled well can be removed.
Third, data-driven approaches such as Principal Component
Analysis (PCA), Wiener filtering and matched filtering and, more
recently, Independent Components Analysis (ICA) have been used
for denoising and artifact rejection (de Weerd, 1981a,b; de Weerd
and Kap, 1981; Lutkenhoner et al., 1984; Makeig, 2002; Makeig etal., 1997; Vigario et al., 2000; Jung et al., 2000; Ikeda and Toyama,
2000; Tang et al., 2002; Sander et al., 2002; James and Gibson,
2003; Lee et al., 2003; Joyce et al., 2004).These approaches haveprovided some benefits but require a subjective choice of many
parameters such as dimensionality, thresholds and component

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

selection. Most of these methods usually do not work well for low
SNR and are computationally prohibitive for multichannel systems
with large number of channels. Furthermore, most of thesealgorithms also provide no principled mechanism for model orderselection. Hence, the selection of the number of spatially distinct
evoked sources and the number of sources of interest must be
based on ad-hoc methods or by expert analysis.
In this paper, we present a novel and powerful approach for the
suppression of interference signals and the separation of signals
from individual evoked sources. This approach is formulated in theframework of probabilistic graphical models with latent variables,
which has been developed and studied in the fields of machinelearning and statistics (Jordan, 1998). In the graphical modelingframework, observed data are modeled in terms of a set of latentvariables, which are signals that are not directly observable. The
dependence of the data on the latent variables is specified by aparameterized probability distribution. The latent variables aremodeled by their own probability distribution. The combined
distributions define a probabilistic model for the observed data.
The model parameters are inferred from data using an expectation
maximization (EM) type algorithm, which is a standard techniquefor performing maximum likelihood in latent variable models.
Problems such as interference suppression and source separation
then translate to the problem of probabilistic (Bayesian) inferenceof appropriate latent variables.
In our case, stimulus-evoked MEG and EEG data are modeledusing a new probabilistic graphical model, which is based on thewell known factor analysis model (Rubin and Thayer, 1982), andon its extension, independent factor analysis (Attias, 1999). The
new model is termed stimulus-evoked independent factor analysis(SEIFA). This model describes the observed data in terms of two
sets of independent latent variables, termed factors. The factors in
the first set represent evoked sources, and the factors in the second
set represent interference sources. The sensor data are generated bylinearly combining the factors in the two sets using two mixing
matrices followed by adding sensor noise. The mixing matrices andthe precision matrix of the sensor noise constitute the SEIFA model
parameters. The interference mixing matrix (the one applied to theinterference factors) and the noise precision are inferred from prestimulus data, whereas the evoked mixing matrix is inferred fromthe post-stimulus data.
The paper is organized as follows. The SEIFA probabilistic
graphical is defined in mathematical terms in SEIFA probabilistic
graphical model. Inferring the SEIFA model from data: a VB-EM
algorithm presents a VB-EM (a generalization of standard
expectation maximization) algorithm (Neal and Hinton, 1998;
Attias, 1999, 2000) for inferring this model from data. Estimating
clean-evoked responses and their correlation matrices provides anestimator for the clean-evoked response, i.e., the contribution of theevoked factors alone to the sensor data, using the model to remove
the contribution of the interference sources. It also breaks thisestimator into the separate contributions from each independent
evoked factor. Moreover, this section presents an automaticallyregularized estimator of the correlation matrix of the clean evoked
response, as well as the correlation matrix of the separate
contribution from each evoked factor. Localization of cleaned
evoked responses demonstrates, using real and simulated data, thatthe algorithm provides interference-robust estimates of the timecourse of the stimulus evoked response. Furthermore, it shows thatusing the regularized evoked covariance in an existing source
localization method improves the performance of that method.

401

Results concludes with a discussion of our results and ofextensions to the SEIFA framework.

SEIFA probabilistic graphical modelThis section presents the SEIFA probabilistic graphical model,
which is the focus of this paper. The SEIFA model describes
observed MEG and EEG sensor data in terms of three types ofunderlying unobserved signals: (1) signals arising from stimulusevoked sources, (2) signals arising from interference sources and (2)sensor noise signals. The model is inferred from data by an algorithm
presented in the next section. Following inference, the model is used
to separate the evoked source signals from those of the interference
sources and from sensor noise, thus providing a clean version of theevoked response. The model further separates the evoked response
into statistically independent factors. In addition, it produces aregularized correlation matrix of the clean-evoked response and ofeach independent factors, which facilitates localization. A graphicalrepresentation of the model is shown in Fig. 1.
Let y in denote the signal recorded by sensor i = 1 : K at timen = 1 : N. We assume that these signals arise from L evoked factors
and M interference factors that are combined linearly. Let x jn
denote the signal of evoked factor j = 1 : L, and let u jn denote thesignal of interference factor j = 1 : M, both at time n. We use theterm factor rather than source for a reason explained below. Let Aijdenote the evoked mixing matrix, and let Bij denote theinterference mixing matrix. Those matrices contain the coefficientsof the linear combination of the factors that produces the data.
They are analogous to the factor loading matrix in the factor
analysis model. Let v in denote the noise signal on sensor i.
Mathematically,
L

M

j¼1

j¼1

yin ¼ ~ Aij xjn þ ~ Bij ujn þ vin

ð1Þ

We use an evoked stimulus paradigm, where a stimulus ispresented at a specific time, termed the stimulus onset time, and isabsent beforehand. The stimulus onset time is defined as n = N 0 + 1.
The period preceding the onset n = 1 : N 0 is termed pre-stimulus
period, and the period following the onset n = N 0 + 1 : N is termed

Fig. 1. SEIFA graphical model. Left column shows the graphical model
for pre-stimulus data. y is the observed spatiotemporal data, B is theinterference factor loading matrix and u the unobserved interference
factors. Right column shows the graphical model for post-stimulus data.
Note the addition of stimulus-evoked factors x, which are modeled as a
mixture of Gaussians defined by states s, and the stimulus-evoked mixing
matrix A. The dashed box encloses all the variables and the nodes
outside the box indicate parameters, all of which are estimated from the
data.

402

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

post-stimulus period. We assume the evoked factors are active only
post-stimulus and satisfy x jn = 0 before its onset. Hence
Bun þ vn ;n ¼ 1 : N0
ð2Þyn ¼
Axn þ Bun þ vn ; n ¼ N0 þ 1 : N
To turn Eq. (2) into a probabilistic model, each signal must be
modeled by a probability distribution. Here, each evoked factor ismodeled by a mixture of Gaussian (MOG) distributions. For factor
j, we have a MOG model with S j components, also termed states,
Sj


  X
N xjn jlj;sj ; mj;sj pj;sj
p xjn ¼

ð4Þ

j¼1

There are two reasons for using MOG distributions, rather thanGaussians, to describe the evoked factors. First, evoked brain
sources are often characterized by spikes or by modulated
harmonic functions, leading to non-Gaussian distributions. Second,previous work on ICA has shown that independent Gaussian
sources that are linearly mixed cannot be separated. Since we aim
to separate the evoked response into contributions from individual
factors, we must therefore use independent non-Gaussian factor
distributions. Third, as is well known, a MOG model with asuitably chosen number of states can describe arbitrary distributions at the desired level of accuracy.
For interference signals and sensor noise, we employ a Gaussianmodel. Each interference factor is modeled by a zero-meanGaussian distribution with unit precision, p ujn ¼ N ujn j0; 1 .
SEIFA describes the factors as independent,pðun Þ ¼

MY 
p ujn ¼ N ðun j0; I Þ

ð5Þ

j¼1

The Gaussian model implies that we exploit only second order
statistics of the interference signals. This contrasts with the evokedsignals, whose MOG model facilitates exploiting higher order statistics, leading to more accurate reconstruction and separation.
The sensor noise is modeled by a zero-mean Gaussian
distribution with a diagonal precision matrix E,pðvn Þ ¼ N ðvn j0; EÞð6ÞFrom Eq. (2) we obtain p( y n |x n ,u n ) = p(v n ) where we substitutev n = y n  Ax n  Bu n with x n = 0 for n = 1 : N 0. Hence, we obtain
the distribution of the sensor signals conditioned on the evoked andinterference factors,N ðyn jBun ; EÞ;n ¼ 1 : N0pðyn jxn ; un ; A; BÞ ¼
N ðyn jAxn þ Bun ; EÞ; n ¼ N0 þ 1 : N
ð7Þ

1

n¼1

A Gaussian distribution over a random vector xhwith mean l and precii
L 1=2
sion matrix  is defined by N ð xjl; LÞ ¼ j 2p
j exp  12 ð x  lÞT Lðx  lÞ .
The precision matrix is defined as the inverse of the covariance matrix.

NYpð xÞ ¼pðxn Þ

n ¼ N0 þ 1

ð3Þ

State s of source j is a Gaussian1 with mean l j,s j and precisionm j,s j , and its probability is p j,s j . We model the factors as mutuallystatistically independent, hence
L
Y
 
p xjn

N

pð yjx; u; A; BÞ ¼ k pðyn jxn ; un ; A; BÞpðuÞ ¼

sj ¼ 1pðxn Þ ¼

SEIFA also makes an i.i.d. assumption, meaning the signals at
different time points are independent. HenceN
Ypðun Þ

ð8Þ

n¼1

where y, x, u denote collectively the signals y n , x n , u n at all time
points. The i.i.d. assumption is made for simplicity and implies thatthe algorithm presented below can exploit the spatial statistics ofthe data but not their temporal statistics.
To complete the definition of SEIFA, we must specify priordistributions over the model parameters. For the noise precision
matrix E, we choose a flat prior, p(E) = const. For the mixing
matrices A, B we use a conjugate prior. A prior distribution istermed conjugate w.r.t. model when its functional form is identicalto that of the posterior distribution (see the discussion below Eq.
(58)). We choose a prior where all matrix elements are independentzero-mean GaussianspðAÞ ¼

Y
N Aij j0; Ei aj

ijpðBÞ ¼YN Bij j0; Ei hj

ð9Þ

ij

and the precision of the ijth matrix element is proportional to thenoise precision Ei on sensor i. It is the E dependence which makes
this prior conjugate. (It can be shown that, in the limit of zero
sensor noise, E Y V, the impact of the prior on the posterior meanof A, B would vanish in the absence of this dependence, which
would be undesirable.) The proportionality constants aj and hj
constitute the parameters of the prior, also known as hyperparameters. Eqs. (8) and (9) together with Eqs. (4), (5) and (7) fully
define the SEIFA model.

Inferring the SEIFA model from data: a VB-EM algorithm
This section presents an algorithm that infers the SEIFA modelfrom data. SEIFA is a probabilistic model with latent variables since
the evoked and interference factors are not directly observable. We
use an extended version of the expectation maximization (EM)algorithm to infer the model from data. This version is termedvariational Bayesian EM (VB-EM) (Attias, 2000).
Standard EM computes the posterior most likely parametervalue given the observed data, aka the maximum a posteriori(MAP) estimate. In contrast, VB-EM considers all possible
parameters values and computes the posterior probability of eachvalue conditioned on the observed data. VB-EM therefore treats
latent variables and parameters on equal footing by computing
posterior distributions for both quantities. One may, however,choose to compute a posterior only over a subset of model
parameters, while computing just a MAP estimate for the rest.
VB-EM is an iterative algorithm, where each iteration consistsof an E-step and an M-step. The E-step computes the sufficient
statistics (SS) of the latent variables, and the M-step computes theS.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

SS of the parameters. (SS of an unobserved variable are quantities
that define its posterior distribution.) The algorithm is iterated toconvergence, which is guaranteed.
The VB-EM algorithm has several advantages compared tostandard EM. It is more robust to overfitting, which can be asignificant problem when working with high-dimensional but
relatively short time series, as we do in this paper. It produces
automatically regularized estimators, such as for the evoked
response correlation matrix, where standard EM produces poorly
conditioned ones. In addition, the variance of the posterior
distribution it computes (essentially the estimator’s variance, or
squared error) provides a measure of the range of parameter values
compatible with the data.
We now describe the VB-EM algorithm for the SEIFA model.
A full derivation is provided in Appendix A.

403

where Ā, B̄, W are computed in the M-step by Eqs. (16) – (18). As
explained in Appendix A, W/Ei is the posterior covariance of row iof AV, l rV, rrV are given in Eqs. (38) and (53).
rr
The covariances Gxx
and Guuof the evoked and interferencerfactors, and their cross-covariance Gxu, conditioned on collective
state r, are obtained by appropriately dividing Gr into quadrants

Gr ¼

Gr
 rxxT
Gxu

Grxu
Gruu


ð14Þ

r
r
where Gxx
is the top left L  L block of Gr, Gxu
is the toprright L  M block, and Guu
is the bottom right M  M block.
These covariances are used in the M-step.
Finally, for a given collective state r, compute its posterior
probability by

E-step
The E-step of VB-EM computes the SS for the latent variables
conditioned on the data. For the pre-stimulus period, n = 1 : N 0, thelatent variables are the interference factors u n . Compute theirposterior mean ū n and covariance A by

ð10Þ

where B̄ are W BB are computed in the M-step by Eqs. ((16) – (18)).
B̄ is the posterior mean of the interference mixing matrix, and W BB
is related to its posterior covariance (specifically, the posteriorcovariance of the ith row of B is W BB /Ei ; see Appendix A).
For the post-stimulus period, n = N 0 + 1 : N, the latent variables
include the evoked and interference factors x n , u n . They also
include the collective state s n of the evoked factors, defined by the
L-dimensional vector s n = (s 1n , s 2n , . . . s Ln ), where s jn = 1 : S j is
the state of evoked factor j at time n. The total number of collective
states is S = C j S j .
To simplify the notation, we combine the evoked andinterference factors into a single vector, and their mixing matrices
into a single matrix. Let LV = L + M be the combined number of
evoked and interference factors. Let AV denote the K  LV matrix
containing A and B, and let x nV denote the LV  1 vector containingx n and u n ,
 
xn
xVn ¼
; AV ¼ ðA BÞ
ð11Þun
The SS are computed as follows. At time n, let r run over all the
S collective states. For a given collective state r, compute theposterior means x̄ rn and ū rn of the evoked and interference factors,
and their posterior covariance Gr, conditioned on r, byx̄rnV ¼ Gr ĀVT Eyn þ rVr lVrÞ
1
Gr ¼ ĀVT E ĀVþ rVr þ KW

ð12ÞHere, as in Eq. (11), we have combined the posterior means ofV , and the posterior means of thethe factors into a single vector x̄ rn
mixing matrices into a single matrix ĀV,x̄rnV ¼




x̄ rn
; ĀV¼ Ā B̄
ū rn



1 pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1
1
1
1
pr jrr jjGr jexp  yTn Eyn þ lTr rr lr  x̄V
rn Gr x̄V
rn
zn
2
2
2
ð15Þ

where z n is a normalization constant ensuring ~r p̄ rn = 1 and the
MOG parameters l r, rr are given in Eq. (38).

Tū n ¼ AB̄ Eyn T
1
A ¼ B̄ EB̄ þ I þ KWBB

p̄rn ¼

ð13Þ

M-stepThe M-step of VB-EM computes the SS for the model
parameters conditioned on the data. We divide the parameters into
two sets. The first set includes the mixing matrices A, B, for which
we compute full posterior distributions. The second set includes thenoise precision E and the hyperparameters matrices a, h, for which
we compute MAP estimates.
Compute the posterior means of the mixing matrices byĀ ¼ Ryx W

B̄ ¼ Ryu W

ð16Þ

where

W¼

Rxx þ aRTxu

Rxu
Ruu þ h1
ð17ÞThe quantities R yx , R yu , R xx , R xu and R uu are posterior
correlations between the factors and the data and among thefactors themselves and are computed below. The hyperparameters
aj , hj are diagonal entries of diagonal matrices a, h.
The covariances W AA and W BB corresponding to the evoked
and interference mixing matrix (see Appendix A), and W AB
corresponding to their cross-covariance, are obtained by appropriately dividing W into quadrants

W¼

WAA
WTAB

WAB
WBB


ð18Þ

where W AA is the top left L  L block of W, W AB is the top rightL  M block, and W BB is the bottom right M  M block.

404

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

Next, use those covariances to update the hyperparameter
matrices A, B bya1 ¼ diagh1 ¼ diag





1 T
Ā E Ā þ WAA
K
1 T
B̄ E B̄ þ WBB
K



ð19Þ

In this section, we present two sets of estimators computed by
the SEIFA model after inferring it from data. The first set of
estimators computes the clean-evoked response and decomposes it
into independent factors. The second set of estimators computes
well-conditioned correlation matrices for the signals obtained by
the first set.
Let z inj denote the individual contribution from evoked factor j
to sensor signal i, and let z in denote the combined contribution
from all evoked factors to sensor signal i. Then

ð20Þ

zinj ¼ Aij xjn

and to update the noise precision matrix E by
E1 ¼



1
T
diag Ryy  ĀRTyx  B̄R yu
N

Posterior means and correlations of the factors

SXp̄ rn x̄ rn

S
X

j¼1

N
X

~n by n u nT ynx̄ Tn

Cj ¼

Let R xx = ~n bx n x nT , R xu = ~n bx n u nT  and R uu = ~n bu n u nT 
denote the evoked – evoked, evoked – interference and interference – interference posterior correlations. ThenC¼

ynū Tn

N
X

S
X

S
X

n¼1

S
X



p̄ rn ū rnū Trn þ Gruu

N
X
n¼1

yn yTnNX T
znj znjNX

zn zTnð28Þ

n ¼ N0 þ 1

ð29Þ

ð23Þ

n ¼ N0 þ 1 r ¼ 1

T

C ¼ ĀRxx Ā þ E1 TrðRxx WAA Þusing the factors covariances (14).
Finally, let R yy denote the data – data correlation
Ryy ¼

ð27Þ

where ā j denotes the jth column of Ā and is a K  1 vector, andthe posterior correlations matrices are given by Eqs. (55), (56) and(60). For the latter correlation matrix, we similarly obtainn ¼ N0 þ 1 r ¼ 1
N
X

z̄inj

j¼1

h  
i
T
j
C̄ ¼ ā j ā j þ E1 ðWAA Þjj ðRxx Þjjp̄ rn x̄ rnū Trn þ Grxu

N0
X
 T

ū nū n þ A þ

LX

Let C̄ j and C̄ denote the estimators of C j and C. This means, as
above, that C̄ j = bC j  and C̄ = bC. For the former, we obtainp̄ rn x̄ rn x̄ Trn þ Grxx

n ¼ N0 þ 1 r ¼ 1
NX

Āij x̄jn ¼

n ¼ N0 þ 1

n¼1

Ruu ¼

ð26Þ

Next, consider the correlation matrix of the evoked response,
which is a required input for localization algorithms such as
beamforming. Let C j denote the correlation of the individual
contribution from evoked factor j, and let C denote the correlationof the combined contribution from all evoked factors. Then

ð22Þ

Rxu ¼

ð25Þ

j¼1

j¼1

n ¼ N0 þ 1
N
X

L
X

z̄in ¼

Let R yx ==denote the data –evoked and data – interference posterior correlations. Then

Rxx ¼

zinj

For the combined contribution we obtain~n by n x nT  and R yu

Ryu ¼

L
X

Let z̄ inj and z̄ in denote the estimators of z inj and z in . This meansthat z̄ inj = bz inj and z̄ in = bz in  where the average is w.r.t. theposterior over A, x. Computing these estimates amounts toobtaining a clean version of the individual contribution from eachfactor and of their combined contribution and removing contributions from interference factors and sensor noise. For the individual
contributions, we obtain

ð21Þ

p̄ rnū rn

r¼1

Ryx ¼

Aij xjn ¼

z̄ inj ¼ Āij x̄jn

r¼1

ū n ¼LX

zin ¼

Here, we compute the posterior correlations, used above,
between the factors and the data and among the factors themselves.
Let x̄ n = bx n  and ū n = bu n  denote the posterior mean of theevoked and interference factors. During the pre-stimulus periodn = 1 : N 0, x̄ n = 0 and ū n is given by Eq. (10). During thepost-stimulus period n = N 0 + 1 : N, they are given byx̄ n ¼

Estimating clean-evoked responses and their correlation
matrices

ð24Þ

ð30Þ

We point out two facts about C̄ j and C̄. First, unlike the signal
estimates which satisfy z̄ in = ~j z̄ inj, their correlation matrices satisfy~j C j m C. Second, and most importantly, all estimated correlation
matrices are well-conditioned due to the diagonal W AA terms.
Hence, the VB-EM approach automatically produces regularized

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

405

correlation matrices. Notice that the correlation matrix obtained
directly from the signal estimates, ~n z̄n z̄nT, is poorly conditioned.

vector corresponding to the minimum eigenvalue of the matrix:[FT (r)Rz1F(r)] (Sekihara et al., 2004).

Localization of cleaned evoked responses

Results

The localization of cleaned sensor responses from the individual evoked factors z nj and from the combined contribution of all
evoked factors z n can be achieved with many algorithms, such as
maximum-likelihood dipole estimation or tomographic reconstruction methods. In this paper, we use adaptive spatial filters that use
data correlation matrices as inputs for localization because thesemethods have been shown to have superior spatial resolution andzero localization bias (Sekihara et al., 2005). Here, we use anadaptive spatial filtering algorithm called the minimum-variancevector beamformer which we briefly review here (Van Veen et al.,
1997; Sekihara et al., 2001).
Let z in denote the magnetic field of the cleaned evoked
response (either from a single factor or from the sum of all evokedfactors) measured by the ith sensor at time n, and let the column
vector z n = [z 1n , z 2n , . . ., z Kn ]T denote the set of cleaned evoked
response data across sensor array. Let Rz denote the correlation
matrix of the cleaned sensor data. For SEIFA, we computecorrelation matrix Rz of this evoked response as described in theprevious section for each evoked factor (C̄ j ) or the summed
contribution from all evoked factors (C̄). Many denoisingalgorithms such as SVD (Sekihara et al., 2002) and JADE
(Cardoso, 1999) do not provide regularized source correlation
matrices. In such cases, one must apply regularization methods
such as Tikhonov or its variants to the matrix ~n z n z nT.
We assume a current dipole at each location whose position,orientation and magnitude are fixed over time. Let the spatial
location of each brain source be denoted by a 3-D vector r = (x,y, z). The dipole moment magnitude of each source at time n isdenoted by q n . The orientation of the source g is denoted by theangles between its moment vector and the x, y and z axes as g x ,g y and g z , respectively. The forward field vector for the xcomponent of a source at r is denoted by f x (r) = [ f 1x (r), f 2x (r),  , f Kx (r)]T, where f kx (r) expresses the kth sensor output inducedby the unit-magnitude source that is located at r and directed in
the x direction. The forward field matrix F(r) = [ f x (r), f y (r),f z (r)] represents the sensitivity of the sensor array at r in all
directions. The forward field vector which represents the sensitivityof the sensor array in the g direction at r is denoted by f (r, g) =F(r)g T.
An adaptive spatial filter estimate of the moment magnitude, q nlocated at r and directed in the g direction uses the following linear
spatial filter operation:

Simulated dataq n ðrÞ ¼ wT ðr; g; Rz Þzn

ð31Þ

The column vector w(r, g, Rz ) represents a set of weights thatcharacterizes the property of the beamformer. In this paper, we use
a linearly constrained minimum variance vector beamformer with
optimized orientation estimates, where the weight vector w(r) for
each voxel is:wðr; g; Rz Þ ¼

R1
z FðrÞgTT
g F ðrÞR1z FðrÞg

ð32Þ

The optimum orientation g opt for each voxel is computed by
maximizing bqnqnT  with respect to g and is given by the eigen-

Computer simulations are shown for a small number of evoked
and background sources, each of which has a sinusoidal timecourse of random frequency (chosen uniformly over a finite range).
The finite length of each evoked source is enforced by modulating
the corresponding sinusoid by a Hanning window having random
length and random placement within the post-stimulus window. In
contrast, the background sources were not modulated by a Hanningwindow. A spherical volume conductor model is assumed, andsimulated sources are assumed to lie in the x = 0 (coronal) plane. A
275 axial gradiometer whole-head MEG system is simulated bypassing the source activity through an estimated forward field
matrix using a random physical location and orientation for eachneural signal and then adding white Gaussian noise to each sensor.
The power of the evoked sources relative to the power of the
background sources (computed in sensor space and averaged overall channels) is referred to as the signal-to-interference ratio (SIR),
which is defined by
1
0 K N
XX
2
ðAi xn Þ C
BC
Bi ¼ 1 n ¼ 1
CðdBÞSIR ¼ 10log10 B
ð33Þ
C
BX
K X
N
@
2A
ðB i u n Þi¼1n¼1

where A i and B i are the i th rows of A and B, respectively.
Likewise, the signal-to-noise ratio (SNR) defined by,
1
0 K N
XX
2
ðAi xn Þ C
BC
Bi ¼ 1 n ¼ 1
CðdBÞ
SNR ¼ 10log10 B
ð34Þ
C
B X
K X
N
A
@
2
vin
i¼1n¼1is used to quantify the power of the evoked sources relative to thepower of the additive sensor noise. A larger SNR is used for shorter
data lengths to emulate the fact that averaging data over trials has
the effect of increasing the stimulus-locked signal power relative tothe power of the sensor noise at the cost of reducing N, the numberof data samples. The SNR for the simulations is fixed at 10 dB
for N = 1000 and at 5 dB for N = 10,000. The stimulus onset,
N 0, is 375 for N = 1000 and is 3750 for N = 10,000. The results
below are given for several different values of SIR.
The comparisons include the proposed method (SEIFA), an
ICA method (JADE), results obtained by selecting the mostenergetic singular-value components (SVD) and the results
obtained when no denoising procedure is used (raw data). For
the simulated results, it is assumed that the number of evoked
sources and background sources is known. Consequently, thenumber of components extracted using SEIFA is explicitly set to L.
ICA is performed using JADE following an SVD truncation, where
the number of components kept during the truncation is L + M.
The denoised sensor signals are based on the L ICA componentsthat have the highest ratio of post-stimulus energy to pre-stimulus
energy. The SVD results are based on keeping the largest L
components. For SVD and JADE, a regularized covariance matrix

406

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

was computed on the denoised sensor data output and used for
beamforming.
Fig. 2 shows results for two configurations of evoked andinterference sources. The left two columns show a simulation
with two evoked sources, L = 2, three interference sources, M =3, with N = 10,000 and SIR = 0 dB and SNR = 5 dB. The righttwo columns show a simulation with L = 2, M = 20, N = 1000,
SNR = 5 dB and SIR = 10 dB. The time axis is represented interms of milliseconds, where the sample frequency is 1 kHz. The
time axis is also shifted so that 0 ms corresponds to the stimulus
onset, N 0 + 1. The first and third panel of the first row in thisfigure show the true locations of the evoked sources, each ofwhich is denoted by , and the true locations of the backgroundsources, each of which is denoted by . The second and fourthcolumns of the first row show the time courses of the evoked
sources as they appear at the sensors. The time courses of theactual sensor signals, which also include the effects of background sources and sensor noise, are shown in the second andfourth columns of the last row. Rows 2 – 4 show the localizationand the time courses of the denoised outputs of the competing&

methods. For M = 3, background sources SEIFA and JADE
provide good localization of the two evoked sources, whereas theraw data and SVD activations have spread to include one or morebackground sources. Furthermore, the raw data include an
additional peak corresponding solely to a background source.
When M is increased to 20, the results for JADE, SVD and theraw data include one or more spurious peaks, and the activation
for JADE becomes much more diffuse. The results for SEIFA,
however, are largely unaffected by the increase in the number ofbackground sources.
Fig. 3 shows the separation results for the single dataset that
corresponds to the left two columns of Fig. 2. Notice that SEIFA
separates the two sources quite well and provides an accurate
localization of each individual evoked source. JADE, which is able
to separate the evoked sources from the background sources as
indicated by Fig. 2, is unable to completely separate the two
evoked sources from each other. The error in the separation is
sufficient so that the localization of each Findependent_ component
produces erroneous results. Notice that, due to the nonlinearity(due to beamforming) involved in constructing the localization

Fig. 2. Simulation results. The left two columns show a simulation with two evoked sources, L = 2, three interference sources, M = 3, with N = 10,000 and
SIR = 0 dB and SNR = 5 dB. The right two columns show a simulation with L = 2, M = 20, N = 1000, SNR = 5 dB and SIR = 10 dB. For all rows, first
and third panels from left show the localization of the stimulus-evoked sources. Second and fourth panel from left show the time courses of the evoked
sources as they appear at the sensors. Top row shows the true location of the evoked sources which is denoted byand the true locations of the
background sources, each of which is denoted by . Second row—results from SEIFA. Third row—corresponding results from JADE. Fourth row—corresponding results from SVD. Fifth row—results from raw data. For M = 3 background sources (left two columns), both SEIFA and JADE provide good
localization of the two evoked factors, whereas the raw data and SVD activations have spread to include one or more background sources. When M is
increased to 20, the results for JADE, SVD, and the Raw Data include one or more spurious peaks and the activation for JADE becomes much morediffuse. Results for SEIFA, however, are largely unaffected by the increase in the number of background sources.

&

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

407

Fig. 3. Separation of individual (estimated) evoked factors for the simulated dataset that corresponds to the left two columns of Fig. 1, for L = 2, M = 3, K =275, N = 10,000. SEIFA separates the two sources quite well and provides an accurate localization of each individual evoked source. JADE, which is able toseparate the evoked sources from the background sources as indicated by Fig. 1, is unable to completely separate the two evoked sources from each other.

maps, the sum of the two individual localizations in Fig. 3 is notguaranteed to equal the single localization of the combined result
shown in Fig. 2. Results are not shown for SVD because of poor
separation performance.

Figs. 2 and 3 correspond to two specific examples. Figs. 4 and
5, on the other hand, show the mean results across 50 Monte Carlotrials as a function of SIR for L = 2 and for four combinations of M
and N; namely, M/N = 3/1000, 20/1000, 3/10,000 and 20/10,000.

Fig. 4. Simulation results using 50 Monte Carlo trials for L = 2, M = 3, K = 275, N = 1000 in the left column and for L = 2, M = 3, K = 275, N = 10,000 in theright column. An SNR of 10 dB is assumed for all simulations, and performance is plotted as a function of the input SIR. For each condition, the location andfrequency of the evoked sources and the interference sources were randomized. Top row—mean output SNIR for SEIFA, JADE, SVD and Raw Data. Middle
row—mean SSNIR for SEIFA and JADE. Bottom row—mean localization error for SEIFA, JADE, SVD and Raw Data.

408

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

Fig. 5. Legend same as for Fig. 3 but for larger number of interference sources M = 20.

For these figures, the amount of denoising is quantified using theoutput signal-to-(noise + interference) ratio (SNIR), which isdefined by10NXð y4in Þ2 C
B X
C
B1 K
n ¼ N0 þ 1
C ðdBÞ
Output SNIR ¼ 10log10 B
C
BK
N
X
@ i¼1
2A
ð yV

y4
Þ
in
in
n ¼ N0 þ 1

ð35Þ
where y *in is the portion of the ith sensor at time n due only to theset of evoked sources and y inV represents the denoised signal at theith sensor found by applying one of the methods above to theobserved (raw) data. Larger values of the output SNIR represent
improved performance, and an infinite value represents perfectdenoising.
The separation metric is given by the (separated-signal)-to(noise + interference) ratio (SSNIR),
01
SSNIR ¼B
L
B
1 X
10log10 B
B
L j¼1
@

1
2  2 N 1 N0

N
X

x4jn xjn

C
C
CðdBÞ
C
A

n ¼ N0 þ 1

ð36Þ
where x*jn is the value of the jth evoked factor at time n, x jn is the jthestimated factor, both x*jn and x jn are normalized to have unit
variance, and the ordering of the factor estimates is such that it
produces the best alignment between the evoked factors and thefactor estimates. This last requirement is needed due to afundamental permutation indeterminacy of BSS. The argument of
the absolute value represents the normalized correlation coefficient
between a particular evoked source and the corresponding source
estimate, which takes values between 1 and +1. Similar to the
output SNIR, the separation metric becomes infinite when thesource estimates perfectly match the true evoked sources. The
output SNIR can be thought of as a collective measure of separation
between the set of evoked sources and the set of backgroundsources/sensor noise, whereas the separation metric provides amean measure of separation between a specific evoked source andall other signals (including other evoked sources). The localization
error is the mean distance in centimeters between the true evoked
source locations and the estimated locations. It is defined byL
1 X
Localization Error ¼
jjr4j  rj jj ðcmÞ
ð37Þ
Lj¼1
where r*j is the true location of the j th evoked source in the y – z
plane, r j represents the estimated location, the estimated locations
are determined as the locations of the L largest peaks that result
from the beamforming reconstruction of the set of denoised
observations, and the estimated locations are paired with the true
source locations in the manner that minimizes the mean distance.
Unlike the preceding two metrics, smaller values of the localization
metric indicate improved performance. Error bars, not shown in the
figures, are on the order of 0.4 dB, 0.9 dB and 0.3 cm for the
output SNIR, separation metric and localization error, respectively.
In terms of denoising performance, SEIFA provides a 5 – 10 dB
improvement over JADE for small N and a 2 – 3 dB advantage over
JADE for large N. Both SEIFA and JADE outperform SVD, which
provides no improvement over the raw data for low SIR and asignificant improvement over the raw data for large SIR. In terms

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

of separation performance, the results for SEIFA far outperform
those of JADE for this set of examples. JADE is able to separate
the background sources from the evoked sources (hence gives
good denoising performance), but it is not always able to separate
the evoked sources from each other (the results shown in Fig. 3 arefairly representative). The Infomax algorithm (Bell et al., 1995),
results not shown, also exhibited poor separation performance onthese data, similar to JADE. These ICA algorithms perform poorly
for the separating of these source presumably due to incorrect
assumptions about the source distributions and due to limited data
length (especially in the case of N = 1000). In terms of localization,
notice that using SVD actually hurts performance. JADE hastrouble for the most challenging case, where N is small and M + L
is large. SEIFA, on the other hand, performs quite well across allconditions and consistently outperforms the other three methods.
Note that we model background sources as Gaussians instead ofbimodal probability distributions corresponding to sinusoidal
sources. Yet, performance of SEIFA is not affected by violationof this assumption. Therefore, according to these data, SEIFA
appears to be quite robust to the i.i.d. assumption of the evoked andbackground sources.
Real data
SEIFA is also applied to real data produced from four
experiments. These experiments are chosen to demonstrate theperformance of SEIFA for denoising, artifact rejection and two
examples of separation of multiple evoked sources. MEG data forthese experiments were collected using a 275-channel, whole-head
MEG system (Omega 2000, VSM MedTech Inc., Port Coquitlam,

409

Canada). The MRI data were collected using either a 1.5 T Signa or
3.0 T Signa Excite by General Electric (Milwaukee, WI). For eachexperiment, either an auditory and/or somatosensory stimulus arepresented for a total of N avg number of trials. The Fraw_ magnetic
field is defined as the trial-averaged data, which are found byaligning the trials based on the timing of the stimulus presentationand computing the mean.
For visualization of the results of these experiments, three kindsof displays are used. First, the time courses of the cleaned evokedfactors before and after denoising are displayed. Second, a contour
map that shows the polarity and magnitude of the denoised and raw
sensor signals in sensor space is shown. The contour plot of themagnetic field on the sensor array corresponds to the mapping of
three-dimensional sensor surface array to points within a circle.
The contour map typically shows the magnetic field profile at aparticular instant of time relative to the stimulus presentation. Thirdand finally, localization of a particular evoked factor or the sum of
all evoked factors are overlaid on the subjects’ MRI. Threeorthogonal projections—axial, sagittal and coronal MRI slices—that highlight all voxels having 80% or more of the energy found in
the maximum-energy voxel (the energy of each voxel is timeaveraged over the entire post-stimulus period) are shown. All
results are shown in the neurological convention. MRI coregistration is done using the NUTMEG software toolkit (http://
www.mrsc.ucsf.edu/~sarang/nutmeg/).
Denoising averages from small number of trials
Results for the first experiment are shown in Figs. 6 and 7. The
stimulus is a 400 ms duration, 1 kHz tone presented binaurally
once every 1.35 s (T100 ms). Fig. 6 shows a representative set of

Fig. 6. Denoising auditory-evoked responses. Each row shows the averaged auditory-evoked response for N avg = 10, 50, 150, 250 respectively. Columns show
the performance of SEIFA, JADE, SVD and Raw Data. Using SEIFA, it can be seen that an evoked response is observable even for averages with very small
number of trials.

410

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

Fig. 7. Denoising performance for auditory-evoked data. The output SNIR is plotted as a function of the number of trials used in the averages. The inset shows
the signals used as a surrogate reference for ground-truth response. Superior performance of SEIFA can be observed for averages from a small number of trials.

the denoised sensor signals (every tenth channel) from each of four
methods as a function of the number of trials used to compute theaverage. The number of components used for SEIFA is L = 1, and
SVD uses the 10 largest singular vectors, and JADE was performed
after preprocessing with SVD. The results for JADE and SVD

assume L = 5 and M = 5. Subsequently, for reconstruction of
denoised signals with SVD and JADE; specific components were
picked to maximize pre – post-stimulus power. The results from all
four methods produce meaningful time courses for large N avg, butthere is a significant deterioration of the results as N avg is reduced.

Fig. 8. Auditory-evoked responses contaminated by eye blink artifacts. Each factor localization plot in this figure and all subsequent ones consist of five panels.
First, localizations are overlaid on the subjects’ MRI. Three orthogonal projections-coronal (top left), sagittal (top right) and axial (center left) MRI slices areshown. In each MR overlay, all voxels having 80% or more of the energy found in the maximum-energy voxel (the energy of each voxel is time-averaged overthe entire post-stimulus period) are shown in white. The bottom row shows the time course of the cleaned evoked factors after denoising are also displayed. A
contour map (center right) that shows the polarity and magnitude of the sensor signals in sensor space is shown for a particular instant of time as indicated by
the vertical dotted line in the bottom panel. (a) The averaged auditory-evoked responses that is dominated by the eye blink artifact. (b) SVD results for sensors
in the left hemisphere, using 2 largest singular vectors. (c) SVD results for sensors in the right hemisphere, using 2 largest singular vectors.

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

Notice that the M100 and the M150 responses for SEIFA are stillclearly distinguishable for N avg as small as 10.
For real data, the time courses at the sensors generated solely by
the evoked sources, y *in, are not known, hence the previously
defined denoising metric is not directly applicable. In an attempt toquantify the performances of the different methods, a filteredversion of the raw data for N avg = 250 is used in place of y *in.
Filtering was done using a zero-phase band-pass filter with cutoff
frequencies of 1 and 100 Hz. Fig. 7 shows the output SNIR as afunction of N avg and in the inset shows the signals used for y *in. For
this example, the output SNIR for SEIFA is nearly flat for N avg >50. SEIFA also produces the best performance of the four methods
for N avg < 225. As can be seen, the SVD method outperforms theothers as N avg becomes large since SVD becomes an ideal methodin this case, i.e., for large SIR.
Artifact rejection
Results for the second experiment are shown in Figs. 8 and 9.
The stimulus is a binaural presentation of a 400 ms duration, 1 kHz

411

tone every 1 s. Data were collected for N avg = 100 trials. The subject
was also asked to blink both eyes every other tone presentation toensure that, for the purpose of this experiment, the eye blink
artifacts occur synchronously with the stimulus. This unrealistic
scenario is chosen for demonstrative purposes only. SVD and
SEIFA are applied separately to the left-half and right-half sensors.
In both cases, SVD uses the two largest singular vectors and SEIFA
assumes L = 2 since two evoked sources are expected, one due tothe eye blink and one due to the auditory stimulus. The eye blink
may be considered to be Fevoked_ for these data due to theexperimental design. Fig. 8a shows the original trial-averaged
magnetic fields recorded at the sensors. Notice that the subject was
anticipating the regularly occurring tones. It is surmised that the
subject started the eye blink early so that her eyes were shut just as
the tone occurred. The auditory response is completely obscuredby the eye blinks, which have a much larger magnitude. Figs. 8band c show the two results for SVD, one for each subset of sensors.
The contour plots are shown for t = 150 ms since this corresponds
to the expected timing of one of the two largest auditory responses

Fig. 9. Localization of SEIFA estimated factors for eye blink data. Legend for each sub-figure is similar to Fig. 7. (a and b) SEIFA results for first and second
factor from sensors in left hemisphere. (c and d) SEIFA results for the first and second factor from sensors in the right hemisphere. In each hemisphere, a factor
can be localized to the eye blink and another to auditory cortex.

412

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

Fig. 10. Localization of SEIFA estimated evoked sources for auditory – somatosensory data in Experiment #3. Three factors extracted using SEIFA for the
combined auditory/somatosensory stimulus are shown. Legend for each sub-figure is identical to Fig. 7. (a and b) Two factors that localize to somatosensory
cortex correspond to and early and late somatosensory response. (c) A third factor localizes to auditory cortex and has the time course consistent with anauditory-evoked response.

and since the eye blink magnitude at this time is relatively small.
The reconstruction using the left-half sensors shows a response in
A1, as desired. The reconstruction using the right-half sensor,
however, does not show a response in A1.
Figs. 9a and c show the reconstructions of the first of the two
factors using SEIFA for the left and right hemispheres, respectively. These are clearly due to the eye blink artifact. Figs. 9b and dshow the second factors for both hemispheres, which both show
activity in A1. Notice that the extracted auditory responses have asignificantly smaller magnitude than the eye blink. This demonstrates the ability of SEIFA to remove artifacts, although it shouldbe mentioned that there are other approaches that could be used forthis particular dataset. Namely, one could use SVD and remove themost energetic singular-value component, which is expected to berelated to the eye blink.

Separation of evoked sources I
To highlight SEIFA’s ability to separately localize evoked factors
and sources, we conducted an experiment involving simultaneouspresentation of auditory and somatosensory stimuli. We expected
simultaneous activation of contralateral auditory and somatosensory
cortex. A pure tone (400 ms duration, 1 kHz, 5 ms ramp up/down)
was presented binaurally with a delay of 50 ms following apneumatic tap on the left index finger. Averaging is performed overN avg = 100 trials triggered on the onset of the tap. Results for thisexperiment are shown in Fig. 10. Results are based on using the righthemisphere channels above contralateral somatosensory and auditory
cortices. The figure shows localization and time course of each of thethree factors extracted by SEIFA. The first two factors localize toprimary somatosensory cortex (SI), however, with differential
latencies. The first factor shows a peak response at a latency of 50

Fig. 11. Localization of SEIFA estimated auditory-evoked sources in Patient F. Legend for each sub-figure is the same as in Fig. 8. (a) SVD beamformer results
using two largest singular vectors using sensors in both hemispheres localize to the center of the head. (b) Localization of the first factor estimated by SEIFA
corresponds to auditory cortex in the left hemisphere and has an early response (85 ms). (c) Localization of the second factor estimated by SEIFA corresponds
to auditory cortex in the right hemisphere with a later response (100 ms).

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

ms, whereas the second factor shows the response at a later latency.
Interestingly, the third factor localizes to auditory cortex, and theextracted time course corresponds well to an auditory evoked
response that is well-separated from the somatosensory response.
Separation of evoked sources II
The patient involved in this example represents an interesting
clinical case in that the auditory response in the right hemisphere isnoticeably delayed relative to the response in the left hemisphere.
This is particularly pronounced for the case of 500 Hz tones, such as
used in this example. The cause of this anomaly is unknown but
may be related to the presence of a brain tumor on one side.
Auditory responses are otherwise normal. Responses to such data
are normally analyzed by beamformer localization one hemisphereat a time. Beamforming of the entire sensor array produces(incorrect) localizations in the center of the brain due to thepresence of highly correlated source activity in both hemispheres.
SEIFA was run on data from the entire sensor array, and results areshown in Fig. 11. Fig. 11a shows the SVD reconstruction, which isbased on the two largest singular vectors. The reconstruction shows
center of the head localization. This is expected since SVD is notable to separate spatially distinct sources, and the activation issymmetric with respect to the center of the head. SEIFA was able tosuccessfully separate the sources in each hemisphere. Figs. 11b – cshow the first and second factors extracted by SEIFA, respectively.
The first factor shows significant activity near the left A1, and thesecond factor shows activity near the right A1. The time course of
these factors clearly shows the early response in the left hemisphere(¨85 ms) when compared with the right (¨100 ms). Hence, SEIFA
is able to separate constituent sources in this experiment.

Discussion
We have presented a novel algorithm for denoising, separationand localization of stimulus-evoked brain activity obtained from
MEG and EEG measurements. Although the results presented here
focus on MEG data, the framework presented here is applicable toEEG data and combined measurements as well. We used a latentvariable probabilistic modeling approach to model sensor measurements as the sum of unobserved non-Gaussian stimulus-evoked
sources, Gaussian interference sources and Gaussian sensor noise.
We used an VB-EM algorithm to infer the latent variables andmodel parameters. We further exploited the pre-stimulus data toinfer parameters of interference sources and sensor noise. We
demonstrate the superior performance of our algorithm with otherbenchmarks both for simulations with high noise and interference
and small data lengths and with real data.
Whereas the framework presented in this paper bears somesimilarities to the Independent Components Analysis (ICA)
approach, there are several important differences. ICA algorithms
learn unmixing matrices and are limited to square noiseless mixingof independent brain sources. In the case of measurement systems
with 275 channels, for good performance, dimensionality reduction
is necessary. Often, SVD methods are used in conjunction withICA algorithms for dimensionality reduction. Signal sub-spacedimensionality estimation can be difficult for low SIR and SNR
data, and dimension reduction will result in erroneous results. In
contrast, SEIFA learns unknown mixing matrices for cases that are
not limited to square-mixing and offers a natural form ofdimensionality reduction with no loss of information. Furthermore,413

it is optimal in non-zero noise and also can exploit pre-stimulus
data to learn the interference model. As we have shown, SEIFA
significantly outperforms JADE in our simulations, especially forshort data lengths and high noise conditions.
One advantage of the algorithm presented here is that it offers aprincipled way of model order selection and to learn evoked source
distributions. In our simulations and data analysis presented so far,
we have assumed that the number of factors L is known. In ouralgorithm, one can use the MAP estimates of the hyperparametersof the mixing matrices to estimate the number of factors by
thresholding. Alternatively, one can determine the number offactors in the data by computing the evidence or marginal
likelihood for different number of factors and then choose L suchthat the marginal likelihood of the data is maximized. This involves
computing Eq. (46) for different model orders and choosing themodel order that maximizes F . We also assume that the source
distributions are known and correspond to damped and windowed
sinusoids. Using the EM algorithm, we can also in principle learn
the means and variances of these mixture of Gaussian distribution
parameters from the data. Model order selection and learning ofsource distributions within the SEIFA framework presented in thispaper are currently under investigation.
SEIFA improves localization of denoised and separated sources
with beamforming by computing regularized data covariancematrices of evoked sources and the noise. Localization of datadenoised using SEIFA can be achieved by algorithms other thanbeamforming. Since SEIFA computes regularized evoked sourceand noise correlation matrices, it can also be used in conjunction
with maximum-likelihood dipole estimation procedures and withMUSIC (Mosher et al., 1992). Such investigations are also
currently underway.
Acknowledgments
This work was funded by NIH grant DC004855-01A1. This
work was supported by NIH Grant R01-DC004855 to SSN.
Appendix A. The VB-EM algorithm
This section outlines the derivation of the VB-EM algorithm
that infers the SEIFA model from data.
A.1. Model
We start by rewriting the model in a form that explicitly includes
the collective states s n of the evoked factors as latent variables. Let r
be an L-dimensional vector denoting a particular collective state: r =(r 1, r 2, . . ., r L ) denotes the configuration where evoked factor j is instate r j = 1 : S j ; there are S = C j S j configurations. Let the Ldimensional vector l r and diagonal L  L matrix rr denote themean and precision, respectively, corresponding to collective stater, and let p r denote its probability. Hence
0

0
1
l1;r1
r1;r1
Bl C
B
B 2;r2 CB
lr ¼ BC; r ¼ B@ N A r @lL;rLpr ¼ p1;r1 p2;r2    pL;rL

1
C
C
C;
A

r2;r2
N
rL;rL

ð38Þ

414

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

We now havepðxn jsn ¼ rÞ ¼ N ðxn jlr ; mr Þ; pðsn ¼ rÞ ¼ pr

ð39Þand from the i.i.d. assumptionNYpð xjsÞ ¼

and since the marginal likelihood p( y) is independent of q,
maximizing F w.r.t. q is equivalent to minimizing the KL distance.
Furthermore, F is upper bounded by log p( y) because the KL
distance is always nonnegative. Hence, any algorithm that successively maximizes F , such as VB-EM, is guaranteed to converge.
A.3. Derivation of VB-EMpðxn jsn Þ

n ¼ N0 þ 1

NYpðsÞ ¼ð40Þ

psnn ¼ N0 þ 1

The full joint distribution of the SEIFA model is given bypð y; x; s; u; A; BÞ ¼ pð yjx; u; A; BÞpð xjsÞpðsÞpðuÞpðAÞpðBÞ

ð41Þ

BF
¼ blog pð y; x; s; u; A; BÞ2  log qð x; s; ujyÞ þ C 1 ¼ 0
Bqð x; s; ujyÞtogether with Eqs. (5), (7), (8), (39) and (40).
A.2. Variational Bayesian inference
The Bayesian approach, as discussed above, treats latentvariables and parameters on equal footing: both are unobserved
quantities for which posterior distributions must be computed. A
direct application of Bayes’ rule to the SEIFA model would
compute the joint posterior over the latent variables x, s, u andparameters A, B
pð x; s; u; A; BjyÞ ¼

1
pð y; x; s; u; A; BÞpð yÞ

VB-EM is derived by alternately maximizing F w.r.t. the twocomponents of the posterior q. In the E-step, one maximizes w.r.t.
the posterior over latent variables q(x, s, u|y), keeping the second
posterior fixed. In the M-step, one maximizes w.r.t. the posterior
over parameters q(A, B |y), keeping the first posterior fixed. Whenperforming maximization, normalization of q must be enforced byadding two Lagrange multiplier terms to F in Eq. (45).
Maximization is performed by setting the gradients to zero

ð42Þ

where the normalization constant p( y), termed the marginal
likelihood, is obtained by integrating over all other variables
XZpð yÞ ¼dx du dA dB pð y; x; s; u; A; BÞð43Þs

However, this exact posterior is computationally intractablebecause the integral above cannot be obtained in closed form.
The VB approach approximates this posterior using a variational technique. The idea is to require the approximate posterior tohave a particular factorized form then optimize it by minimizing
the Kullback – Leibler (KL) distance from the factorized form to
the exact posterior. Here, we choose a form which factorizes thelatent variables from the parameters given the data

BF
¼ blog pð y; x; s; u; A; BÞ1  log qðA; BjyÞ þ C2 ¼ 0
BqðA; BjyÞ
ð47Þ
where C 1, C 2 are constants depending only on the data y and noton the variables that are used in the E- and M-steps. bI1 denotes
averaging only w.r.t. q(x, s, u|y), and bI2 denotes averaging only
w.r.t. q(A, B |y). Hence, the posteriors are given byqð x; s; ujyÞ ¼qðA; BjyÞ ¼

1
exp ½ blog pð y; x; s; u; A; BÞ2 
Z1

1
exp ½ blog pð y; x; s; u; A; BÞ1 
Z2

where Z 1, Z 2 are normalization constants.
A.4. E-step
It follows from Eq. (48) that the posterior over u, x, s factorizes
over time and has different pre- and post-stimulus forms,
qðu; x; sjyÞ ¼

N0
Yqðun jyn Þ I

n¼1

pð x; s; u; A; BjyÞ , qð x; s; u; A; BjyÞ ¼ qð x; s; ujyÞqðA; BjyÞð44ÞIt is worth emphasizing that (1) beyond the factorization
assumption, we make no further approximation when computing
q and (2) the factorized form still allows correlations among x, s, u,
as well as among the matrix elements of A, B, conditioned on thedata.
Rather than minimize the KL distance directly, it is convenient
to start from an objective function defined by
XZ
F ½q ¼dx du dA dB qð x; s; u; A; BjyÞ½log pð y; x; s; u; A; BÞs log qð x; s; u; A; BjyÞ

ð45Þ

It can be shown that
F ½q ¼ log pð yÞ  KL½qð x; s; u; A; BjyÞjjpð x; s; u; A; BjyÞ

ð46Þ

ð48Þ

N
Y

qðun ; xn ; sn jyn Þ

ð49Þ

n ¼ N0 þ 1

It also follows that in the pre-stimulus, period q(u n |y n ) is
Gaussian in u n , and, in the post-stimulus period, q(u n , x n , s n |y n ) is
Gaussian in u n , x n for a given s n . To see this, consider log q(x, s,u|y) in Eq. (48) and observe that it is a sum over n, where the nthelement depends only on x n , u n and the dependence is quadratic.
For the pre-stimulus period, we obtainqðun jyn Þ ¼ N un jū n ; A1ð50Þ

with mean ū n and covariance matrix A given by Eq. (10). (Onefirst obtains A = (bBT EB + I) 1 and then performs the averageusing Eq. (61).) For the post-stimulus period, we write the posterior
in the formqðun ; xn ; sn jyn Þ ¼ qðxn ; un jsn ; yn Þqðsn jyn Þ

ð51Þ

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

so w ji = AijV . Then, the posterior over each row is Gaussian

The first component on the right hand side is Gaussianqðxn ; un jsn ¼ r; yn Þ ¼ qð xVnjsn ¼ r; yn Þ ¼ N xVnjx̄Vrn; G1
r

415

ð52Þ

K
Y

qðA; BjyÞ ¼ qð AVjyÞ ¼
N wi j w̄ i ; Ei W1

ð58Þ

i¼1

with mean x̄rn
V and covariance matrix Gr1 given by Eq. (12) (as for
A above, one first obtains Gr = (bAVT E AV + I) 1 then applies Eq.
(61).) l rV and rrV are defined as the LV  1 mean vector and LV  LV
diagonal precision matrix, respectively, for collective state r,
obtained from m r and Nr of Eq. (38) by padding with M zeros

lVr ¼

lr
0M




rVr ¼



rr

ð53Þ

0M

The state posterior is given by q(s n = r|y n ) = p̄ rn (15), which is
obtained by some further algebra.
It is useful to make explicit the correlations among the factors
implied by their posteriors ((50), (52)). For the pre-stimulus period,
we obtain
bun uTn  ¼ ū n ū Tn þ Að54Þwith mean w̄ ji = Ā ij computed by Eq. (16). The precision matrix
Ei W  1 is computed using Eq. (17). To see this, consider log q(A, B |y)in Eq. (48) and observe that it is a sum over i, where the ith element
depends only on the ith rows of A, B and the dependence isquadratic.
It is now evident that p(A, B) of Eq. (9) is indeed a conjugate
prior. Rewriting it in the form
pðA; BÞ ¼ pðAVÞ ¼

KYN wi j0; Ei aVÞ

where aV is a diagonal matrix with the hyperparameter matrices A,
h on its diagonal, shows that its functional form is identical to thatof the posterior (58), with W 1 replacing aV.
It is useful to make explicit the correlations among the elementsof the mixing matrices implied by their posterior (58). They are
bAijV AklV  = ĀijV Ā klV + yik W jl /Ei , or in terms of A, B

For the post-stimulus period, we obtain conditioning on thecollective state, bx nV |s n = r = x̄ rn
V and bx nV x nV T |s n = r = x̄ rn
V x̄ rn
VT + G r.
In terms of xn , u nbAij Akl  ¼ Ā ij Ā kl þ dik

bxn jsn ¼ r ¼ x̄ rnbBij Bkl  ¼ B̄ ij B̄ kl þ dik1ðWAA Þjl
Ei

1ðWBB Þjl
Ei
1bAij Bkl  ¼ Ā ij B̄ kl þ dik ðWAB ÞjlEibun jsn ¼ r ¼ ū rn

ð59Þ

i¼1

ð60Þ

where we used Eq. (18). It follows thatbxn xTn jsn ¼ r ¼ x̄ rn x̄ Trn þ Grxx

TbAT EA ¼ Ā EĀ þ KWAA
bun uTn jsn ¼ r ¼ ū rn ū Trn þ Gruu

T

bBT EB ¼ B̄ EB̄ þ KWBB
bxn uTn jsn ¼ r ¼ x̄ rn ū Trn þ Grxu

ð55Þ

where we have used Eqs. (13) and (14).
To obtain the posterior means and correlations ((21), (23)), wemust sum over the collective state since
XZ
X
dxnV qðsn ¼ rÞqð xnV jsn ¼ rÞxnV ¼
p̄ rn bxnV jsn ¼ rbxnV ¼rbxnV xVnT  ¼

XZ

r

X

p̄ rn bxnV xVnT jsn ¼ r

ð56Þ

r

Eqs. (55) and (56) lead to the results (21) and (23).

blog pð yjx; u; A; BÞ þ log pðAÞ þ log pðBÞIt follows from Eq. (48) that the parameter posterior factorizes
over the rows of the mixing matrices and correlates their columns.
Let w i denote a column vector containing the ith row of the
combined mixing matrix AV = (A, B)

and set its derivative w.r.t. L to zero.


B 
B
AV ¼ @



1

w
w2

wK

1


C
C
A


ð57Þ

ð62Þ

where the averaging is w.r.t. the posterior q. Next, compute thederivative of this expression w.r.t. a, h and set it to zero. The
solution of the resulting equation is Eq. (19). It is easier to firstcompute the derivative and then apply the average. Similarly, toobtain the update rule for the noise precision (20), observe that thepart of F that depends on E is

A.5. M-step

0

ð61Þ

which are needed for Eqs. (10) and (12).
To obtain the update rules for the hyperparameters (19), observethat the part of the objective function F (45) that depends on A, B
isblog pðAÞ þ log pðBÞdxnV qðsn ¼ rÞqð xnV jsn ¼ rÞxnV xVnT

r¼

bAVT EAV ¼ ĀVT EĀVþ KW

ð63Þ

References
Attias, H., 1999. Independent factor analysis. Neural Comput. 11 (4),803 – 851.
Attias, H.A., 2000. Variational-Bayesian framework for graphical models.
Adv. Neural Inf. Process. Syst. 12, 209 – 215.
Baillet, S., Mosher, J.C., Leahy, R.M., 2001. Electromagnetic brain
mapping. Signal Process. Mag. 18, 14 – 30.

416

S.S. Nagarajan et al. / NeuroImage 30 (2006) 400 – 416

Baryshnikov, B.V., Van Veen, B.D., Wakai, R.T., 2004. Maximumlikelihood estimation of lowrank signals for multiepoch MEG/EEG
analysis. IEEE Trans. Biomed. Eng. 51 (11), 1981 – 1993.
Bell, A.J., Sejnowski, T.J., 1995. An information maximisation approachto blind separation and blind deconvolution. Neural Comput. 7 (6),1129 – 1159.
Berg, P., Scherg, M., 1991. Dipole modelling of eye activity and its
application to the removal of eye artefacts from the EEG and MEG.
Clin. Phys. Physiol. Meas. 12 (Suppl. A), 49 – 54.
Bijma, F., de Munck, J.C., Huizenga, H.M., Heethaar, R.M., 2003. A
mathematical approach to the temporal stationarity of background noise
in MEG/EEG measurements. NeuroImage 20 (1), 233 – 243.
Bijma, F., de Munck, J.C., Bocker, K.B., Huizenga, H.M., Heethaar, R.M.,
2004. The coupled dipole model: an integrated model for multiple
MEG/EEG data sets. NeuroImage 23 (3), 890 – 904.
Cardoso, J.-F., 1999. High-order contrasts for independent component
analysis. Neural Comput. 11 (1), 157 – 192 (January).
de Munck, J.C., Bijma, F., Gaura, P., Sieluzycki, C.A., Branco, M.I.,
Heethaar, R.M., 2004. A maximum-likelihood estimator for trial-to-trial
variations in noisy MEG/EEG data sets. IEEE Trans. Biomed. Eng. 51(12), 2123 – 2128.
de Weerd, J.P., 1981. Facts and fancies about a posteriori wiener filtering.
IEEE Trans. Biomed. Eng. 28 (3), 252 – 257.
de Weerd, J.P., 1981. A posteriori time-varying filtering of averaged evoked
potentials: i. introduction and conceptual basis. Biol. Cybern. 41 (3),211 – 222.
de Weerd, J.P., Kap, J.I., 1981. A posteriori time-varying filtering of
averaged evoked potentials: ii. mathematical and computational aspects.
Biol. Cybern. 41 (3), 223 – 234.
Ebersole, J.S., 1997. Magnetoencephalography/magnetic source imaging
in the assessment of patients with epilepsy. Epilepsia 38 (Suppl. 4),
S1 – S5.
Hmlinen, M., Hari, R., Ilmoniemi, 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.
Ikeda, S., Toyama, K., 2000. Independent component analysis for noisy
data—MEG data analysis. Neural Netw. 13 (10), 1063 – 1074.
James, C.J., Gibson, O.J., 2003. Temporally constrained ICA: anapplication to artifact rejection in electromagnetic brain signal analysis.
IEEE Trans. Biomed. Eng. 50, 1108 – 1116.
Jordan, M.I. (Ed.), 1998. Learning in Graphical Models. The MIT Press,
Cambridge, MA.
Joyce, C.A., Gorodnitsky, I.F., Kutas, M., 2004. Automatic removal of eye
movement and blink artifacts from EEG data using blind component
separation. Psychophysiology 41, 313 – 325.
Jung, T.P., Makeig, S., Humphries, C., Lee, T.W., McKeown, M.J., Iragui,V., Sejnowski, T.J., 2000. Removing electroencephalographic artifacts
by blind source separation. Psychophysiology 37 (2), 163 – 178.

Lee, P.L., Wu, Y.T., Chen, L.F., Chen, Y.S., Cheng, C.M., Yeh, T.C., Ho,
L.T., Chang, M.S., Hsieh, J.C., 2003. Ica-based spatiotemporalapproach for single-trial analysis of postmovement MEG beta synchronization. NeuroImage 20 (4), 2010 – 2030.
Lutkenhoner, B., Wickesberg, R.E., Hoke, M., 1984. A posteriori wiener
filtering of evoked potentials: comparison of three empirical filter
estimators. Rev. Laryngol. Otol. Rhinol. (Bord.) 105 (2), 153 – 156.
Makeig, S., 2002. Response: event-related brain dynamics—Unifying brain
electrophysiology. Trends Neurosci. 25 (8), 390.
Makeig, S., Jung, T.P., Bell, A.J., Ghahremani, D., Sejnowski, T.J., 1997.
Blind separation of auditory event-related brain responses into independent components. Proc. Natl. Acad. Sci. U. S. A. 94 (20), 10979 – 10984.
Mosher, J.C., Lewis, P.S., Leahy, R.M., 1992. Multiple dipole modeling andlocalization from spatio-temporal MEG data. IEEE Trans. Biomed. Eng.
39 (6), 541 – 557.
Neal, R.M., Hinton, G., 1998. A view of the EM algorithm that justifies incremental, sparse and other variants. In: Jordan, M.I. (Ed.), Learning in Graphical Models. Kluwer Academic Publishers, Dordrecht,
pp. 355 – 368.
Rubin, D.B., Thayer, D.T., 1982. EM algorithms for ML factor analysis.
Psychometrika, 69 – 76.
Sander, T.H., Wübbeler, G., Lueschow, A., Curio, G., Trahms, L., 2002.
Cardiac artifact subspace identification and elimination in cognitive
MEG data using timedelayed decorrelation. IEEE Trans. Biomed. Eng.
49, 345 – 354.
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(7), 760 – 771.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., Miyashita, Y.,
2002. Application of an meg eigenspace beamformer to reconstructing
spatio-temporal activities of neural sources. Hum. Brain Mapp. 15,
199 – 215.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., 2004. Asymptotic
snr of scalar and vector minimum-variance beamformers for neuromagnetic source reconstruction. IEEE Trans. Biomed. Eng. 51 (10),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 (4), 1056 – 1067.
Tang, A.C., Pearlmutter, B.A., Malaszenko, N.A., Phung, D.B., Reeb, B.C.,
2002. Independent components of magnetoencephalography: localization. Neural Comput. 14 (8), 1827 – 1858.
Van Veen, B.D., van Drongelen, W., Yuchtman, M., Suzuki, A., 1997.
Localization of brain electrical activity via linearly constrained minimum
variance spatial filtering. IEEE Trans. Biomed. Eng. 44 (9), 867 – 880.
Vigario, R., Sarela, J., Jousmaki, V., Hamalainen, M., Oja, E., 2000.
Independent component approach to the analysis of EEG and MEG
recordings. IEEE Trans. Biomed. Eng. 47 (5), 589 – 593.