www.elsevier.com/locate/ynimg
NeuroImage 37 (2007) 102 – 115

A probabilistic algorithm integrating source localization and noise
suppression for MEG and EEG data
Johanna M. Zumer, a,b Hagai T. Attias, c Kensuke Sekihara, d and Srikantan S. Nagarajan a,b,⁎aBiomagnetic Imaging Lab., Department of Radiology, University of California, San Francisco, San Francisco, CA 94143-0628, USA
UCSF/UC Berkeley Joint Graduate Group in Bioengineering, University of California, San Francisco, San Francisco, CA 94143-0628, USAc
Golden Metallic Inc., PO Box 475608, San Francisco, CA 94147, USAd
Department of Systems Design and Engineering, Tokyo Metropolitan University, Tokyo 191-0065, Japanb

Received 31 December 2006; revised 29 March 2007; accepted 19 April 2007
Available online 13 May 2007
We have developed a novel probabilistic model that estimates neuralsource activity measured by MEG and EEG data while suppressing theeffect of interference and noise sources. The model estimates
contributions to sensor data from evoked sources, interference sources
and sensor noise using Bayesian methods and by exploiting knowledge
about their timing and spatial covariance properties. Full posterior
distributions are computed rather than just the MAP estimates. In
simulation, the algorithm can accurately localize and estimate the timecourses of several simultaneously active dipoles, with rotating or fixedorientation, at noise levels typical for averaged MEG data. The
algorithm even performs reasonably at noise levels typical of an
average of just a few trials. The algorithm is superior to beamforming
techniques, which we show to be an approximation to our graphicalmodel, in estimation of temporally correlated sources. Success of thisalgorithm using MEG data for localizing bilateral auditory cortex,
low-SNR somatosensory activations, and for localizing an epileptic
spike source are also demonstrated.
© 2007 Elsevier Inc. All rights reserved.
Keywords: Magnetoencephalography; Electroencephalography; Inverse
methods; Bayesian inference; Denoising

Introduction
Mapping the spatiotemporal neural activity of the entire brain is
an important problem in basic neuroscience research. It is also
clinically important for patients with brain tumors and epilepsy,
both in localizing regions important for cognitive function and for
identifying epileptogenic brain regions. Such brain mapping⁎ Corresponding author. Biomagnetic Imaging Lab., Department ofRadiology, University of California, San Francisco, San Francisco, CA
94143-0628, USA. Fax: +1 415 502 4302.
E-mail address: sri@radiology.ucsf.edu (S.S. Nagarajan).
URL: http://bil.ucsf.edu (S.S. Nagarajan).
Available online on ScienceDirect (www.sciencedirect.com).
1053-8119/$ - see front matter © 2007 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2007.04.054

procedures are useful to guide neurosurgical planning, navigation,and resection.
Many noninvasive techniques have emerged for functional
brain mapping, such as functional magnetic resonance imaging(fMRI) and magnetoencephalography (MEG). Although fMRI is
the most popular method for functional brain imaging with highspatial resolution, it suffers from poor temporal resolution since itmeasures blood oxygenation level-dependent (BOLD) signals withfluctuations in the order of seconds. These BOLD signals are also
indirect measures of neural activity which might not accuratelyreflect neural activity, especially in regions of altered vasculature.
However, dynamic neuronal activity has fluctuations in thesubmillisecond time scale that can only be directly measured withmagnetoencephalography (MEG) and/or electroencephalography(EEG). MEG data are the measurements of tiny magnetic fields
surrounding the head while EEG data are the measurements ofvoltage potentials using an electrode array placed on the scalp.
The past decade has shown rapid development of whole-head
MEG/EEG sensor arrays and of algorithms for reconstruction ofbrain source activity from MEG and EEG data, termed source
localization. All existing methods for brain source localization arehampered by the many sources of noise present in MEG/EEG data.
The magnitude of the stimulus-evoked neural sources is on theorder of noise on a single trial, and so typically 50–200 averaged
trials are needed in order to clearly distinguish the sources above
noise. This limits the type of cognitive questions that can be asked
and is prohibitive for examining processes such as learning that canoccur over just one or several trials. Needing to average trials istime consuming and therefore difficult for a subject or patient tohold still or pay attention through the duration of the experiment.
Gaussian thermal noise or Gaussian electrical noise is present at the
MEG or EEG sensors themselves. Background room interferencesuch as from powerlines and electronic equipment can beproblematic. Biological noise such as heartbeat, eyeblink or othermuscle artifact can also be present. Ongoing brain activity itself,
including the drowsy-state alpha (∼10 Hz) rhythm can drown out
evoked brain sources. Finally, many localization algorithms have

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

difficulty in separating neural sources of interest that have
temporally overlapping activity.
Noise in MEG and EEG data is typically reduced by avariety of preprocessing algorithms before being used by source
localization algorithms. Simple forms of preprocessing include
filtering out frequency bands not containing a brain signal ofinterest. Additionally and more recently, ICA algorithms havebeen used to remove artefactual components, such as eyeblinks(Jung et al., 2000). More sophisticated techniques have also
recently been developed using graphical models for preprocessing prior to source localization (Nagarajan et al., 2005, 2006).
Therefore, current algorithms for source localization from MEG
and EEG data typically use a two-stage procedure — the firstfor noise/interference removal and the second for source
localization.
This paper presents a probabilistic modeling framework for
MEG/EEG source localization that estimates Source Activity using
Knowledge of Event Timing and Independence from Noise andInterference (SAKETINI). The framework uses a probabilistic
hidden variable model that describes the observed sensor data interms of activity from unobserved brain and interference sources.
The unobserved source activities and model parameters areinferred from the data by a variational Bayesian Expectation
Maximization algorithm. The algorithm then creates a spatiotemporal image of brain activity by scanning the brain, inferring themodel parameters and variables from sensor data, and using themto compute the likelihood of a dipole at each grid location in thebrain.
We first describe the generative model for the data. We
complete specification of the model for post-stimulus sources by
including prior distributions and compute the unknown quantities
learned from the data. We then describe the model for learninginterference and noise sources from pre-stimulus data. We finish
the Methods by showing that an established source localization
method, the minimum variance adaptive beamformer (MVAB)(Sekihara et al., 2001), is an approximation of our framework. In
the Results section, we show performance of SAKETINI relative to
MVAB and sLORETA (Pascual-Marqui, 2002) both in localization
ability and in time course estimation for MEG data. We show theeffect of number of sensors and time points for all three methods.
We further show the proposed method's performance applied to areal auditory-evoked MEG dataset, a low-SNR somatosensory
MEG dataset, and an epileptic spike MEG dataset. We conclude
with a discussion on the model order and inputs to the algorithm,
extensions of the model, and its relationship to other methods inthe literature. A preliminary report of this work was shown inZumer et al. (2007).
Methods
Probabilistic model integrating source localization and noise
suppression
This section describes the generative model for the data. We
assume that the MEG/EEG data have been collected such that
stimulus onset or some other experimental marker indicated the“zero” time point. Ongoing brain activity, biological noise,background environmental noise, and sensor noise are present inboth pre-stimulus and post-stimulus periods; however, the evokedneural sources of interest are only present in the post-stimulus time
period. We therefore assume that the sensor data can be described103

as coming from four types of sources: (I) evoked source at aparticular voxel (grid point), (II) all other evoked sources not atthat voxel, (III) all background noise sources with spatial
covariance at the sensors (including brain, biological, or environmental sources), and (IV) sensor noise. We first infer the model
describing source types (III) and (IV) from the pre-stimulus data,then fix certain quantities (described below in next subsection) andinfer the full model describing the remaining source types (I) and(II) from the post-stimulus data (described below in the third
subsection). After inference of the model, a map of the source
activity is created, as well as a map of the likelihood of activity
across voxels.
Let yn denote the K × 1 vector of data from sensor data for time
point n, where K is the number of sensors (typically ∼ 200). Timeranges from − Npre:0:Npost − 1 where Npre(Npost) indicates thenumber of time samples in the pre-(post-)stimulus period. The
generative model for data yn is
yn ¼

Bun þ vnr
F r srn þ Aqr xq
n þ Bun þ vn

n ¼ Npre ; N ;  1
n ¼ 0; N ; Npost  1

ð1Þ

The K × q forward lead field matrix F r represents the physical (andlinear) relationship between a dipole source at voxel r for each
dipole orientation q along a coordinate basis and its influence on
sensor k = 1:K (Sarvas, 1987). The lead field F r is calculated fromknowing the geometry of the source location to the sensor
location, as well as assuming the type of conducting medium inwhich the source lies. The human head is most commonly
approximated as a single-shell sphere volume conductor for MEG
data; for EEG data, a more detailed model taking into account thetissue conductance boundaries is usually needed. In the mostgeneral case, q = 3 for all three possible directions of coordinate
bases of a source dipole. In the case of the single-shell sphere, theradial component of source dipoles contributes nothing to MEG
sensors, thus q = 2. If there is knowledge of subject-specificcortical anatomy, the source may be constrained to be perpendicular to the gray matter surface, thus q = 1. Throughout the rest of
this paper, we use the single-shell model with q = 2 for both
simulations and real data from MEG.
The source activity srn is a q × 1 vector of dipole strength attime n for the voxel r. The K × L matrix Ar and the L × 1 vector
xrn represent the post-stimulus mixing matrix and evoked nonlocalized factors, respectively, corresponding to source type (II)discussed above, where the \r superscript indicates for all voxels
not at voxel r. We leave these superscripts off in the rest of thispaper for clarity. The K × M matrix B and the M × 1 vector un
represent the background mixing matrix and background factors,
respectively, corresponding to source type (III) discussed above.
The K × 1 vector vn represents the sensor-level noise, which isassumed to be drawn from a Gaussian distribution with zeromean and precision (inverse covariance) defined by the diagonal
matrix λ. All quantities depend on r in the post-stimulus period
except for B, un and λ, which will be learned from the prestimulus data, and fixed as the other quantities are learned foreach voxel. Note, however, that the posterior update for un(represented by ūn) does depend on the voxel r. The graphicalmodel is shown in Fig. 1. This generative model becomes aprobabilistic model when we specify prior distributions, as
described in the next two sections (see Jordan (1999) and Bealet al. (2003) for additional information on related graphicalmodels).

104

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

Fig. 1. (Left) SAKETINI graphical model. Variables dependent on time are inside dotted box; parameters independent of time outside dotted box. Values incircles are unknown and learned from the model, and values in squares are known. After applying a simpler model to the pre-stimulus data to learn the sensor
noise precision and interference mixing matrix, these quantities are then known and placed in squares. (Right) Representation of factors influencing the data
recorded at the sensors. In orange, a post-stimulus source at the voxel of interest, focused on by the lead field F. In red, other post-stimulus sources not at thatparticular voxel. In green, all background sources, including on-going brain activity, eyeblinks, heartbeat, and electrical noise. In blue, thermal noise present ineach sensor.

Separation of background sources learned from pre-stimulus data
We learn the background mixing matrix B and sensor noise
precision λ = diag(λ1,…,λK) from the pre-stimulus data using avariational Bayesian factor analysis model. All variables, parameters and hyperparameters are hidden and are learned from thepre-stimulus data. We assume Gaussian prior distributions on the
background factors and sensor noise; we assume a flat prior on thesensor noise precision. We further assume that the backgroundfactors are independent and identically distributed (i.i.d.) acrosstime. All normal distributions specified in this paper are defined bytheir mean and precision (inverse covariance). The following isassumed for n = −Npre:− 1:pðuÞ ¼ j pðun Þ;pðun Þ ¼ N ðun j0;IÞ;pðvÞ ¼ j pðvn Þ;pðvn Þ ¼ N ðvn j0;kÞ;nnpðyju; B;kÞ ¼ j pðyn jun ; B; kÞ;nð2ÞpðkÞ ¼ const:

ð3Þpðyn jun ; B;kÞ ¼ N ðyn jBun ; kÞð4Þ

We define the notation pðun Þ ¼ N ðun jA;DÞ to mean that theprobability distribution of un is a normal Gaussian distribution withmean μ and precision (inverse covariance) Δ. To complete thespecification of this model, we need prior distributions on themodel parameters. We use a conjugate prior for the backgroundmixing matrix B, as follows,pðBÞ ¼ j pðBkm Þ;
kmpðBkm Þ ¼ N ðBkm j0;kk bm Þð5Þ

where βm is a hyperparameter over the mth column of B and λk is
the precision of the kth sensor. The matrix β = diag(β1, … ,βM)
provides a robust mechanism for automatic model order selection,
so that the optimal size of B is inferred from the data through β (forgeneral reference, see Neal (1996)).

In general, for any “free” probability density function (pdf) q,
the “free energy” F (so-called in statistical physics) is written asZ
F ðq;HÞ ¼

du dB qðu;Bjy;HÞ½log pðy;u;BjHÞ  logqðu;Bjy;HÞ¼ Eqðu;Bjy;HÞ ½lðy;u;BjHÞ þ HqRpðu;Bjy;HÞ
þ logpð yjHÞ
¼ du dB qðu;Bjy;HÞ log
qðu;Bjy;HÞ
¼ log pðyjHÞ  KL½qðu;Bjy;HÞjjpðu;Bjy;HÞ¼ LðHÞ  KL½qðu;Bjy;HÞjjpðu;Bjy;HÞ

ð6Þ
ð6Þ

where Θ is the set of all model parameters, in this case Θ = {β, λ},and Hq is the entropy of q(u,B|y,Θ). The Kullback–Leibler
divergence between distributions q and p, defined KL(q||p) =∫dx q(x)[log q(x) − log p(x)], is always non-negative and can equal
zero if the two distributions are equal; thus, F is a lower bound to
the log likelihood function L. Exact inference on this model is
intractable using the true joint posterior over the backgroundfactors u and mixing matrix B. We desire the marginal conditional
posterior estimates of the background factors and mixing matrix,
given the data. We could choose to estimate just the maximum aposteriori (MAP) estimates, but this ignores the full distribution
including its precision. We could also choose to estimate points ofthe distribution extensively, using, for example, Markov Chain
Monte Carlo methods (Jun et al., 2005; Gelman and Rubin, 1996);
however this is computationally costly as well as dependent on thesampled points. Instead, we choose to factorize the distribution
assuming conditional independence of the factors and parameters,also termed the mean field approximation.
pðu;BjyÞcqðu;BjyÞ ¼ qðujyÞqðBjyÞ

ð7Þ

This approximation restricts the posterior to a product of factor
distributions but allows the solution to be computed analytically.
F is equal to log p(y|Θ) when the approximation in Eq. (7) istrue, thus making the KL-divergence zero. The variational Bayesian
Expectation Maximization (VB-EM) algorithm (Attias, 1999;
Ghahramani and Beal, 2001; Beal, 2003) iteratively maximizes Fwith respect to (w.r.t.) each factorized distribution to, at least, a local
maximum of F , alternating w.r.t. the posteriors q(u|y) and q(B|y). InJ.M. Zumer et al. / NeuroImage 37 (2007) 102–115

the E-step, F is maximized w.r.t. q(u|y), keeping q(B|y) constant,and the sufficient statistics of the hidden variables are computed. In
the M-step, F is maximized w.r.t. q(B|y), keeping q(u|y) constant,and the MAP estimate of the parameters and hyperparameters is
computed. The q(u,B|y) = q(u|y)q(B|y) that maximizes F is the
distribution that minimizes the KL-divergence, thus making themean field approximation as close to equivalence as possible under
the constraint implied by Eq. (7). A derivation of a similar model is
described in the Appendix of Nagarajan et al. (2006) and thevariational approximation has been shown to be a powerful andflexible tool (Jordan et al., 1999). It can further be shown that, given
our previous assumptions of Gaussian prior distributions, the formof the posterior q(u|y) must also be Gaussian, by showing that thefunction log q is quadratic in u.
In the E-step, maximizing F w.r.t. q(u|y) holding q(B|y) fixed
gives the following posterior estimate update for the factors.
qðujyÞ ¼ jqðun jyn Þ;1nTū n ¼ g B̄ kyn ;qðun jyn Þ ¼ N ðun jū n ;gÞ
T

g ¼ B̄ kB̄ þ Kw

1

ð8Þ

þI

For all variables in this paper, a bar over the variable indicates its
posterior mean. In the M-step, we compute the full posterior
distribution of the background mixing matrix B, including its
precision matrix ψ, and the MAP estimates of the noise precision λand the hyperparameter β.
qðBjyÞ¼N ðBjB̄;wÞ
B̄ ¼ Ryu w; 


1 T
1
B̄ kB̄ þ w ;b ¼ diagK

w ¼ ðRuu þ bÞ11diag Ryy  B̄RTyu
k1 ¼
N

N
X

yn yTn ; Ryu ¼

n¼1

N
X

yn ū Tn

ð10Þ

n¼1

In the stimulus-evoked paradigm, the source strength at eachvoxel is learned from the post-stimulus data. The backgroundmixing matrix B and sensor noise precision λ are fixed, after
having been learned from the pre-stimulus data, as described in theprevious section. We assume that those quantities remain constantthrough the post-stimulus period and are independent of source
location. We assume Gaussian prior distributions on the source
factors and interference factors. We further make the assumptionthat both of these factors are i.i.d. across time.
The source factors have prior precision given by the q × q matrix
Φ, which relates to the strength of the dipole in each of q directions.
We have used the superscript r here as a reminder that s is for only
one voxel at a time, although in general we will leave off thissuperscript.
npðsrn Þ ¼ N ðsrn j0; Ur Þ

ð11Þ

The interference and background factors are assumed to haveidentity precision (uninformed prior).
pðxÞ ¼ j pðxn Þ;npðuÞ ¼ j pðun Þ;npðxn Þ ¼ N ðxn j0; IÞ;pðun Þ ¼ N ðun j0; IÞpðAÞ ¼ j pðAkl Þ;pðAkl Þ ¼ N ðAkl j0; kk al Þ

ð13Þ

kl

We now specify the full model:pðyjs; x; u; A; B; kÞ ¼ j pðyn jsn ; xn ; un ; A; B; kÞ;npðyn jsn ; xn ; un ; A; B; kÞ ¼ N ðyn jFsn þ Axn þ Bun ; kÞ

ð14Þ

Exact inference on this model is intractable using the joint posterior
over the evoked non-localized factors x and mixing matrix A; thus avariational Bayesian approximation for the posteriors is again used.
pðs; x; u; AjyÞcqðs; x; u; AjyÞ ¼ qðs; x; ujyÞqðAjyÞ

ð15Þ

All variables, parameters and hyperparameters are unknown and arelearned from the data. We learn the hidden variables and parametersfrom the post-stimulus data, iterating through each voxel across thebrain, using the VB-EM algorithm. Since maximizing the log p(y|Θ), for Θ = {α, Φ} would be mathematically intractable, we
maximize a lower bound to log p(y|Θ) defined by the free energy
F in the following equation.
Z
F ðq; HÞ ¼ dudx ds dA qðs; x; u; Ajy; HÞ½log pðy; s; x; u; AjHÞ log qðs; x; u; Ajy; HÞ¼ log pðyjHÞ
 KL½qðs; x; u; Ajy; HÞjjpðs; x; u; Ajy; HÞ
ð16Þ

Localization of evoked sources learned from post-stimulus datapðsr Þ ¼ j pðsrn Þ;

We also use a conjugate prior for the interference mixing matrix A,
where α = diag(α1,… ,α L) is a hyperparameter, similar to theexpression for the background mixing matrix, which helps in modelorder selection for determining the size of A from the data.

ð9Þ

The matrices, such as Ryu, represent the posterior covariance
between the two subscripts.
Ryy ¼

105

ð12Þ

In the E-step, the posterior distribution of the factors given thedata is computed:
qðxnVjyn Þ ¼ N ðxnVjx̄ nV;CÞ; x̄ n ¼ C1 A¯ VT kyn ;C ¼ A¯ VT k A¯ Vþ KW þ I V

ð17Þ

where we define:
0 1
s̄ n
B C
x̄ nV¼ @ x̄ n A; A¯ V¼ ðF A¯ B̄ Þ;
0

ū n

U
B
I V¼ @ 0010 0CI 0 A;0 I

00B
W ¼ @0

0
WAA

1
0
C
0A

0

0

0

ð18Þ
ð18Þ

In the M-step, we update the posterior distribution of theinterference mixing matrix A including its precision ΨAA. Note thatthe lead field F is fixed and known based on the geometry of thesensors relative to the head (same quantity assumed known by
beamforming and other approaches) and B̄ was learned and fixedfrom the pre-stimulus data. The sensor noise precision λ is also kept
fixed from the pre-stimulus period. The MAP values of thehyperparameter α and source factor precision Φ are learned here
from the post-stimulus data.
A¯ ¼ ðRyx  FRsx  B̄Rux ÞWAA ;1
WAA ¼ ðRxx þ aÞ1 ; U1 ¼ Rss ;N
1 T
A¯ kA¯ þ WAAa1 ¼ diag
K

ð19Þ

106

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

Explicit definitions for all combinations of posterior covariances areomitted, but some are shown here for example:

Rss ¼

N
X

X
s̄ n s̄ Tn þ N ss ; Rsx ¼

n¼1

Xs̄ n x̄ Tn þ N sx

W ¼ X1 F T Y ;
ð20Þ

n¼1

where ∑ = Γ−10P

is specified as:
P

ss

X

N
X

sx

P 1
su

B
C
BXT P P C
C
¼B
xx
xu C
B sx
@
A
XT XT Psu

xu

estimated the model parameters A, B, and λ, the MAP estimate ofthe dipole source time series is s̄ n = Wyn, where

ð21Þ

uu

In each iteration of VB-EM, the marginal likelihood isincreased. The variational likelihood function (the lower boundon the exact marginal likelihood) is given as follows:

X ¼ FT Y F þ U

ð24Þ

The Φ term exists from the prior distribution of sn. It can be shownthat Eq. (24) is equivalent to Eq. (23).
In the infinite data limit, the data covariance satisfies Ryy =FΦ− 1FT + ϒ− 1. Its inverse is found, using the matrix inversion
lemma, to be R−yy1 = ϒ − ϒFΩ− 1FTϒ. Hence, we obtain
1
1 T
T
T
F T R1yy ¼ ðI  F Y FX ÞF Y ¼ UX F Y

ð25Þ

where the last step used the expression for Ω. Next, weapproximate Ω ≈ FTϒF, which is equivalent to the prior term on
sn having zero precision. We then obtain
W cðF T Y FÞ1 F T Y ¼ ðF T Y FÞ1 XU1 d UX1 F T Y
1 T 1¼ ðF T R1yy FÞ F Ryy ¼ WMVAB

ð26Þ

where the last step uses Eq. (25).
Lr ¼

N 
X


N
jkjjU j 1log
yT kyn  x̄ nVrT Cr x̄ nVr2
jCr j
2 n¼1 nK
þ logjar jjWr j2r

Simulations: generation and analysis
ð22ÞThis likelihood function is dependent on the source voxel r, since
we update the variables for each voxel, and thus a map of thelikelihood across the brain can be displayed. Furthermore, animage of the source power estimates and the time course of activity
at each voxel could be plotted.
We note that the computational complexity of the proposed
algorithm is on the order O(KLNS), roughly equivalent to a single
dipole scan, which is on the order O(N(K2+S)). These are muchsmaller than the complexity of a multi-dipole scan which is on theorder O(NSP) where P is the number of dipoles, and if S representsroughly several thousand voxels.
We further note that the number of hidden variables to be
estimated is less than the number of data points observed, thus not
posing problems for estimation accuracy. See Fig. 7 for performance
as a function of data points.
Relationship of SAKETINI to minimum variance adaptive
beamforming
Minimum variance adaptive beamforming (MVAB) is one of
the best performing source localization techniques (Sekihara et al.,2001, 2002, 2005; Küçükaltun-Yildirim et al., 2006). MVAB
estimates the dipole source time series by ŝn = WMVAByn, where1 T 1WMVAB ¼ ðF T R1yy FÞ F Ryy

ð23Þ

Thus, MVAB also has computational complexity equivalent to asingle-dipole scan, on the order O(N (K2+S)). MVAB is known tofail for two sources that are highly correlated in time (Sekihara etal., 2002), although methods exist to correct this issue (Dalal et al.,2006; Brookes et al., 2007). In this section, we derive that MVAB
is an approximation to SAKETINI.
We start by rewriting Eq. (1) as yn = Fsn + zn, where zn is termed
the total noise and is given by zn = Axn + Bun + vn. It has mean zeroand precision matrix ϒ = (AAT + BBT + λ− 1)− 1. Assuming we have

SAKETINI was tested in a variety of scenarios, including many
simulations with real brain noise and real MEG datasets.
Performance of the method is compared to the MVAB andsLORETA (Pascual-Marqui, 2002). Simulations and real data wereanalyzed using NUTMEG (Neurodynamic Utility Toolbox for
MEG) (Dalal et al., 2004), a toolbox developed using MATLAB(MathWorks, Natick, MA, USA), obtainable from http://bil.ucsf.
edu. NUTMEG is useful for coregistration of fiducial points to a
structural MRI, selection of volume-of-interest, computation of
forward field, filtering and other denoising preprocessing methods,
as well as a variety of source reconstruction methods, including
MVAB, sLORETA and now SAKETINI.
Simulations were created using a variety of realistic source
configurations reconstructed on a 5-mm voxel grid. A single-shell
spherical volume conductor model was used to calculate theforward lead field (Sarvas, 1987). Simulated datasets wereconstructed by placing Gaussian-damped sinusoidal time courses
at specific locations inside a voxel grid based on realistic headgeometry. Sources were set to be active only during a post-stimulus
period, which always comprised 62.5% of the total data available,
while the remaining 37.5% was pre-stimulus data. Typically 700
total data points were used, unless specified otherwise.
In one set of simulations, termed “simulated interference” cases,
we created background activity drawn from the statistical distributions assumed by the model. The simulated sources mentioned above
were projected to the sensors using the computed forward field.
Noise drawn from a Gaussian distribution was added to each sensor
across all time points. Furthermore, to simulate “ongoing brain
activity,” time courses were drawn from a Gaussian noise distribution
and were placed in 30 random locations throughout the brain voxel
grid, occurring in both pre- and post-stimulus periods. Their activity
was projected onto the sensors and added to both Gaussian sensor
noise and source activity. These simulated background brain sources
add noise to the sensors in a spatially correlated manner.
In order to test simulation performance using data with more
realistic (and unknown) statistical distributions, we created a second
set of simulations termed “real brain noise.” Real MEG sensor data
were collected from a CTF MEG System with 275 axial gradio-

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

meters while a human subject was alert but not performing tasks or
receiving stimuli. This real background data thus include real sensor
noise plus real ongoing brain activity that could interfere withevoked sources and adds spatial correlation to the sensor data. Since
throughout this paper we work with averaged data, these real data
were binned into 100 trials of 700 data points each and averaged. We
varied the signal to noise ratio (SNR) and the corresponding signal tonoise-plus-interference ratio (SNIR). SNIR is calculated from theratio of the sensor data resulting from sources only to the sensor datafrom noise plus interference.
Simulation methods: average performance
In addition to some illustrative examples, we examine average
results over many realizations of simulations in which certain
source parameters were varied. Two different source configurations were used: one with sources more near the surface as
depicted in Fig. 3 and the other configuration with deeper
sources. The orientation of the source was fixed in half thesimulations and allowed to rotate over time in the other half. The
correlation of two of the three sources with each other was set tobe 0%, 95%, or 100%; the third source was always uncorrelatedwith the other two sources. Either simulated interference or real
brain noise was used as background noise. SNIR was set at 5 dB,
0 dB, or −5 dB. Each combination of parameters was tested for
10 different randomly generated source time courses and source
orientations. Thus, a total of 720 simulations were run using all
combinations of simulation parameters.
Performance was measured in two ways: localization ability
and estimation of time course. To assess localization ability, it isimportant to take into account source strength, source localization
error, and presence of false positives. Thus we used the ROC(receiver-operator characteristic) method modified for brain
imaging results as suggested by Darvas et al. (2004), which is ameasure of hit rate versus false positive rate. The free-response
ROC (FROC) curve in particular allows for multiple “hits” perimage (Bunch et al., 1978).
We define a “local peak” as a voxel that is greater in value thanits 26 three-dimensional neighbors. We define a “hit” as a local
peak that is within a specified distance of the true location andabove a certain threshold. A “miss” is defined as a true source
location that has no “hit” within the specified distance. A “false
positive” is a local peak above a certain threshold but further than
the specified distance from a true source location. A “truenegative” is any voxel that is none of the above.
FROC curves are generated by varying the threshold andallowable distance error,pthus
theand
ﬃﬃﬃ varyingp
ﬃﬃﬃ tradeoff of sensitivity
pﬃﬃﬃ
specificity. We used 5 * 3mm, 10 * 3mm, or 15 * 3mm for theallowable localization error from a true source. We varied thethreshold to be 30%, 50%, 70% or 90% of the maximum value inthe whole image. We thus recorded a hit rate (HR) and false
positive rate (FR) for each of 12 combinations of threshold/error
for each of the 720 simulations.
Since these HR versus FR points do not increase monotonically,as they would if the threshold was the only criterion varied, wechose to use the measure of A′. A′ is a way to approximate the area
under the ROC curve for one HR/FR point (Snodgrass and Corwin,1988). The larger the area under the ROC, the better the method is
performing, since this means a higher HR relative to FR for
specified thresholds/localization errors. For each simulation, thetwelve computed A′ values were averaged to give one A′ value per

107

simulation. The SAKETINI likelihood map was used as the spatial
map to test localization. For MVAB and sLORETA, the power map(sum of squares of post-stimulus time points) was used.
Simulation methods: effects of amount of data
Previous studies have shown advantage of sensor arrays withlarger number of channels (Hamalainen et al., 1993). We sought todetermine how many sensors were required to prevent performancedegradation. Likewise, increased amount of data points usuallyleads to improved estimation of unknown quantities and therefore
we tested performance with less data available.
To test the effect of the number of sensors, simulations werecreated similarly to those discussed above with three uncorrelated
sources. Two values of SNIR were created using real brain noise:
0 dB and −10 dB. Ten different realizations of source time course
and orientation were tested for each case. All simulations discussed
previously were created using the full 275 channel array from the
CTF system. We tested the method's performance when only arandom subset of sensors were selected, using 150, 74, 37, 25 or 15
sensors.
To test the effect of the number of data points, the full set of 275
channels were used but reduced the available amount of data
points. All previous simulations have used 700 total data points,
where 62.5% were in the post-stimulus period. We kept the ratio ofdata points in the post-stimulus period the same but reduced thetotal number to 300, 200, 150, 100, or 50 points.
Collection of real data
Several real datasets were analyzed with the proposed method
and compared to existing methods. For all data, the 275-channel
CTF MEG System in a magnetically shielded room was used tocollect data. All healthy subjects and one patient with epilepsy
gave written, informed consent to participate in each study,
according to UCSF institutional review board approval.
Auditory
The first dataset was obtained by presenting 120 repetitions of a
1-kHz tone binaurally to a healthy subject, at an intertrial intervalof 1.4 s. The trials were averaged locked to stimulus onset. This
auditory stimulus is known to invoke bilateral auditory cortex to be
active simultaneously, causing problems for the MVAB's ability tolocalize the auditory sources.
Somatosensory
We next examine a somatosensory dataset in which the localization of primary somatosensory cortex is relatively easy for all
methods when many trials are available to average. A small
diaphragm was placed on the subject's right index finger and was
driven by compressed air. The stimulus was given 256 times every
500 ms. However, if we limit the available data to only a small
subset of trials, the lower SNR can become limiting for all
methods. We first applied SAKETINI, MVAB and sLORETA tothe average of all 256 trials to assess performance for the standard(high) SNR case. We then applied all three methods to the averageof only the first 5 trials. To further test if the performance wasconsistent across other sets of just 5-trial averages, we applied thethree methods to the 5-trial average of trials 6–10, 11–15, and 16–

108

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

20. We then averaged the results of these four different results. Any
location found consistently will show up in the average.
Epileptiform spikes
We next tested the proposed method on its ability to localizeinterictal spikes obtained from a patient with epilepsy. No sensory
stimuli were presented to this patient in this dataset, which wascollected in the same MEG device described above. A Registered
EEG/Evoked Potential Technologist marked segments of thecontinuously collected dataset which contained spontaneous
spikes, as well as segments that clearly contained no spikes. One
segment of data with a spike marked at 400 ms, as well as three
additional spikes in the 800 ms segment, was used here as the“post-stimulus” period and a separate, spike-free, segment of equal
length was used as the “pre-stimulus” period.
Results
We first report results from an example simulation. Then we
describe performance for the averages of the simulations discussed
above, including varying the amount of sensors or time points
available. We finish by demonstrating performance in the real
datasets discussed above.
Simulation results
Individual example simulation results
We first examine localization error for a single active source
placed randomly within the voxel grid space. Fig. 2 shows
localization error of SAKETINI relative to the MVAB. The largest
peak in the SAKETINI likelihood map was found and the distancefrom this point to the true source was recorded. Each data point isan average of 20 realizations of the source configuration, with error
bars showing standard error. This simulation was performed for avariety of SNIR's and for both cases of simulated and real brain

Fig. 2. Localization error of SAKETINI compared to MVAB. A single
source was placed randomly in a realistic brain volume and its SNIR was
varied. Error was measured as the three-dimensional distance from the
maximal peak in the reconstruction to the true location.

noise. This figure clearly shows that the error in localization issmaller for SAKETINI (black) than for MVAB (green).
We next examine the performance of the proposed model whenthree active sources are present. Fig. 3 shows performance in one
example with three uncorrelated sources. On the left is thelikelihood map and on the right are the estimated time courses(green) of the three sources compared to the true time courses(blue). The likelihood map correctly locates all three sources andthe time courses are estimated well.
We further test performance when two of the three sources arehighly correlated in time. In one example simulation, Fig. 4 shows
the SAKETINI likelihood map on the left compared to the MVAB
power map on the right. The likelihood map detects the presence of
all three sources, whereas the MVAB fails to find the left source
and instead finds a false source in the middle between the twocorrelated sources.
Simulation results: average performance
We now show performance of the proposed methods according
to the metrics of A′ (area under ROC curve) and time course
estimation, as defined in the Methods. Fig. 5(a) plots A′ for eachmethod, for each value of source correlation and SNIR, and forboth types of interference (simulated and real). SAKETINI (blacklines) has a larger average A′ than MVAB (green) and sLORETA(red) for nearly all values of SNIR and source correlation,
indicating overall better localization ability.
The other main test of performance was the ability to estimate
the source time course. The estimated time courses for all threemethods were obtained from the true source locations, regardlessof whether their respective localization maps found that source as a“hit.” The correlation of the true time course with the estimatedtime course was computed for each simulation and the averages areplotted in Fig. 5(b). SAKETINI estimates the time course better
than MVAB for all simulation types and better than sLORETA innearly all simulation types.
Simulation results: effects of amount of data
Fig. 6 shows simulation performance resulting from reduced
number of sensors. Fig. 6(a) shows the A′ metric for localization
performance for SAKETINI relative to MVAB and sLORETA,
where Fig. 6(b) shows estimation accuracy of the source timecourse. The top row in both (a) and (b) is for SNIR = 0 dB and thebottom row is for SNIR = −10 dB. SAKETINI's performance is
not degraded for as few as 74 sensors for both values of SNIRand can even perform well using as few as 25 sensors for SNIR =0 dB.
Fig. 7 tests the proposed method against the other methods withdecreased number of time points available. Fig. 7(a) shows that
SAKETINI has a larger A′ value relative to MVAB and sLORETA
for all amounts of data tested and for both SNIRs. However,
SAKETINI's ability to estimate the source time course (Fig. 7(b))
decreases with less than roughly 150 data points total, even though
MVAB and sLORETA do not show this sensitivity. Like otherminimum-norm techniques, sLORETA is not data-dependent, andthus is not expected to show sensitivity to number of data points.
The MVAB is dependent on the data to provide an estimate of thedata covariance matrix. Since the simulations in both top andbottom rows are with relatively high noise (SNIR = 0 dB and− 10 dB, respectively), the data covariance estimate might notchange much with decreased data, since it is already noisy (notethe time course correlation does not reach above 0.5 for any

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

109

Fig. 3. SAKETINI's performance in simulated example with three uncorrelated sources. (Left) Blue circles and squares indicate location of true sources. Intensityof map corresponds to normalized likelihood map. (Right) Blue lines indicate simulated time series for each of the three source locations. Green lines indicate
proposed method's estimates of the source time series at those three locations.

number of data points tested at SNIR = − 10 dB). However,
SAKETINI relies on an initial estimate of interference sources fromthe pre-stimulus data; perhaps this estimate is not as robust when
only 30 or less pre-stimulus data points are present to estimate the 30
interference sources projected to 275 sensors, allowing for somecross-talk of interference with “evoked” time estimates.
Results from real data
Auditory
The well-characterized auditory-evoked M100 peak was observed and the sensor map for this latency is shown in Fig. 8(a).
Bilateral auditory cortex is clearly activated. As mentioned above,
MVAB performance deteriorates when neural sources are highly
correlated in time and, in general, bilateral auditory cortex is
known to be correlated (Dalal et al., 2006). Fig. 8(b) shows the
MVAB performance on this dataset: the MVAB fails to find either
left or right auditory cortex and instead finds a false source in thecenter. However, both SAKETINI and sLORETA are able to findbilateral auditory cortex (Figs. 8(c) and (d)). Additionally shownin Fig. 8(d) are the time courses of the SAKETINI source
estimates extracted from the peak voxel on each side. The
correlation of the left with the right source estimates was found tobe 0.93.

Somatosensory
The left panel of Fig. 9(a) shows typical somatosensory-evoked
MEG data with the largest peak at 50 ms, expected to be comingfrom primary somatosensory cortex in the posterior wall of thecentral sulcus. The next three panels of Fig. 9(a) show localization
performance of SAKETINI, MVAB and sLORETA. All threemethods accurately localize activity to the contralateral primarysomatosensory cortex. However, performance changes when only 5
trials are used in the average. The left panel of Fig. 9(b) shows thesensor data averaged over trials 1–5 of the same somatosensory
dataset. The next three panels of Fig. 9(b) show errors in localizationin all methods. SAKETINI shows the least error, relative to the peaklocation found using all 256 trials. We note that other averages of 5
trials showed varied performance, but that, when averaging four
different sets of 5-trial averages together, SAKETINI showed
localization closest to primary somatosensory cortex, as shown in
Fig. 9(c), whereas MVAB and sLORETA mislocalize this source.
Epileptiform spikes
Fig. 10 shows the proposed method's performance on estimating
single spikes relative to a spike-free “pre-stim” period. Fig. 10(a)shows the raw sensor data for the segment containing the marked
spike. Fig. 10(b) shows the location of the equivalent-current dipole(ECD) fit to 20 spikes from this patient; the centroid of this cluster ofFig. 4. Example of SAKETINI and MVAB for correlated source simulation. Blue squares indicate source locations that are highly temporally correlated. Bluecircle indicates simulated source location that was not correlated with the other two sources. (Left) SAKETINI likelihood map. (Right) MVAB power map.

110

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

Fig. 5. Performance of SAKETINI relative to MVAB and sLORETA for variety of simulated datasets. Each data point is an average of 40 simulations, consistingof two different source locations and either a fixed or rotating source orientation. Standard errors were less than 0.05 for all points (not shown). (a) A measure ofarea under ROC curve A′ is plotted for the three methods as a function of SNIR for both simulated and real brain interference. See text for discussion of the A′
metric. (b) The correlation of the estimated with the true time course is plotted for each method.

ECD fits would normally be used clinically. Fig. 10(d) shows the
SAKETINI likelihood map based on the data in Fig. 10(a); the peakis in clear agreement with the standard ECD localization. The
maximum voxel of the SAKETINI likelihood map is 6.3 mm fromthe ECD fit of that particular spike and is 3.2 mm from the ECD
cluster centroid. Fig. 10(c) shows the time course estimated for the
likelihood spatial peak. The spike at 400 ms is clearly seen; thiscleaned waveform could be of use to the clinician in analyzing peak
shape (information usually not available from the ECD). Finally,
Figs. 10(e) and (f) show a source time course from a randomly
selected location far from the epileptic spike source (shown withcross-hairs on bottom right plot), in order to show the low noise level

and to show lack of cross-talk onto source estimates elsewhere. In
addition to this single spike analysis with SAKETINI, the sensor
data of all 20 spikes were first averaged, then examined by both
ECD and SAKETINI.
Typically this average is only used clinically if the ECD of single
spikes gives large parameter fit errors. The difference between the
ECD and SAKETINI localizations based on this average spike data
was only 3.1 mm. Both of these also corresponded well to the 20spike ECD cluster centroid mentioned above, with differences of
only 4.7 mm (ECD) and 3.1 mm (SAKETINI).
Furthermore, the cluster of 20 ECD localizations corresponded
well to the electrocorticography mapping performed in this patient.

Fig. 6. (a) A′ and (b) time course estimation as a function of the number of MEG sensors for simulated data with 3 uncorrelated sources. The top row shows
SNIR = 0 dB and bottom row shows SNIR = − 10 dB using real brain noise. Error bars represent standard error.

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

111

Fig. 7. (a) A′ and (b) time course estimation as a function of the number of total data points for simulated data with 3 uncorrelated sources. The top row of eachshows SNIR = 0 dB and bottom row shows SNIR = −10 dB using real brain noise. Error bars represent standard error.

Intraoperative electrodes were spaced 1 cm apart and the cluster
resided within this 1 cm error. This brain region was chosen forresection and this patient was seizure-free after the operation. This
indicates that the brain region localized pre-surgically andnoninvasively by the ECD cluster (and consistent with the
SAKETINI localization) correctly specified the epileptogenic zone.

Discussion
We have described the proposed method in a graphical model
framework, which is a powerful and flexible technique fordescribing probabilistic dependencies between observed andunobserved quantities. We have chosen a scanning-based method

Fig. 8. Performance of methods using a real auditory-evoked MEG dataset. (a) Sensor map showing bilateral activation at M100 latency. (b) MVAB power map,showing failure to localize sources to auditory cortex. (c) sLORETA power map localizes bilateral auditory cortex. (d) SAKETINI likelihood map localizes
activity in bilateral auditory cortex. Time series from left and right peaks are also shown, where intensity is normalized to maximum value across all voxels andtime points.

112

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

Fig. 9. Performance of methods using real somatosensory data as a function of the number of trials. Left column shows sensor data averaged over varied numberof trials, while remaining columns show localization performance of SAKETINI, MVAB and sLORETA. Row (a) shows performance of the three methods
applied to the average of all 256 trials. Row (b) shows the localization performance to the average of only the first 5 trials. In order to show performance overother subsets of 5-trial averages, the spatial maps in row (c) are spatial averages of the localization of 4 different 5-trial averages. See Methods for details.
Crosshairs in localization maps show peak location within “active” voxels at the slice of peak location, where the threshold for “active” was defined at 90% of the
maximum for all maps.

to formulate the problem, rather than to solve the full tomographic
problem which is very ill parameterized. The variational Bayesian
Expectation Maximization algorithm is used to solve for the valuesof the unknown quantities that maximize the marginal log
likelihood. A map of the maximized likelihood is used to depict
localized sources and estimates of their time course are extracted
from those spatial peaks.
We have shown that the proposed model is not as sensitive totemporally correlated sources as the standard formulation of
MVAB. However, it is possible to reduce the MVAB's dependence
on correlated sources through a modified weight matrix computed
subject to additional constraints, if a rough idea of the location ofsources is known (Dalal et al., 2006).
Model order and initialization
As the number of MEG and EEG channels has increased inrecent years, the ability to accurately localize sources throughout
the brain has increased (Vrba et al., 2004). However, performing
calculations of high-dimensional data, such as inverting a data

covariance matrix, becomes more difficult and can lead to errors.
Meanwhile, the dimensionality of the underlying neural activity
remains the same. Thus, many variations of PCA and ICA havebeen used on MEG/EEG data for removal of noise/artefactual
components as well as for data dimension reduction (Jung et al.,
2000; Ikeda and Toyama, 2000). Factor analysis also aims to
reduce the dimensionality of the data to a linear mixture of factorsthat best account for the data while accounting for noise at thesensor level. We have used an extended version of factor analysis
to partition which factors are event-related activity and which
factors are background interference.
All methods which perform dimension reduction need acriterion for choosing the reduction number. Using PCA, a plotof eigenvalues can often give a reasonable intuition for thedimension of “signal” in the data. ICA has no ordering ofcomponents. In the method proposed here, there are two variables
affecting model dimension: number of background interference
factors (u) and number of non-localized evoked factors (x).
The model itself has two built-in ways of assisting in thedimension choice. First, the hyperparameters over the columns of

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

113

Fig. 10. SAKETINI applied to interictal spike data from a patient with epilepsy. (a) Segment of data contains several spikes marked by a technician. No stimulus
was presented to the patient in this data. (b) Cluster of Equivalent Current Dipole (ECD) localizations to 20 spikes are shown from this patient. The middle
column shows SAKETINI time course estimation (c) and likelihood map (d) for this segment of data, showing spatial agreement with ECD. The right column
shows SAKETINI time course (e) from a voxel (f) randomly selected far away from the location of the spike localization, indicating lack of crosstalk of the spikes
into normal tissue. The time course intensities in panels c and e are normalized to maximum value across all voxels and time points.

the mixing matrices A and B are updated in the M-step; thus their
final values can be plotted and used as indication of which columns(and therefore which factors) contribute more, similar to a plot of
eigenvalues for PCA. As long as the number of dimensions chosen
is larger than the true dimension, the extra dimensions can be“zeroed out” by the model as those corresponding hyperparameters
decrease.
A second way to examine model order is by learning the model
for several values of model order and then comparing the final
maximized likelihood value. In general, the greater the modelorder, the easier to explain the data, thus the likelihood will
increase for greater model order. However, there is usually an
inflection point in the plot of likelihood versus model order, again
similar to a plot of eigenvalues for PCA, around which the ideal
model order resides.
All methods which do not have a closed form solution require
initialization of the values to be iteratively updated. We havefound that choice of initialization can change the final results
somewhat but not largely, and so we did not extensively examine
these effects. After finding one method of initialization thatworked well in a few test simulations, that set was used for all
results shown. We note that in testing data with only sensor noise(no simulated or real brain interference at all), we do see a greater
dependence on initialization values. This is worth mentioning
since many source reconstruction algorithms are only tested on
simulated sources plus sensor noise, thus comparison of thismethod with other methods on sensor-noise-only data should only
be done with this in mind. However, since this type of data does
not follow the assumptions of the model that the pre-stimulus

data has spatially covarying noise, and is therefore not realistic,
we chose not to focus the results of this paper on this type ofdata.
A consideration of algorithms that iteratively update values isthe longer computation time. SAKETINI infers the model values
for each voxel and convergence typically occurs after about 10–20
VB-EM iterations, which takes roughly one second per voxel on astandard Linux personal computer with 2.0-GHz processor.
Changing the number of channels or time points also affects
computation time.
Extensions of model
The graphical model proposed here is one variation of a class ofmodels that could be used to perform source reconstruction for
MEG and EEG data. Certainly many aspects of this model can be
modified in order to relax assumptions made here. One assumption
is that the background interference sources are stationary across the
pre-stimulus and post-stimulus periods. If background is nonstationary, the model could be modified to include a term similar toΦ which would influence the posterior update for ū. This would
occur, for instance, if ongoing brain activity decreased as a result of
the stimulus/event onset.
If the background is stationary, but the source activations vary
drastically over the post-stimulus period, it might help to break up
the post-stimulus data into partitions of relatively similar activity.
Specifically, the data could be decomposed into time–frequency
windows where inference of this model on each window would
result in a separate likelihood map.

114

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

In all results shown, we used a single-shell spherical conductor
model for generation of the forward field F. Certainly moresophisticated forward fields could be computed using multisphere
methods or BEM models, which could further improve performance for all methods. Distributed sources could be modeled
explicitly using a forward field computed from patch bases(Limpiti et al., 2006). However, errors in coregistration of MEG
or EEG data to the structural MRI can affect any computation of
forward field. If one wishes to restrict the source activity to one
orientation based on cortical surface, errors in this segmentation
can lead to errors in orientation estimation. Localization of EEG
data has larger errors since tissue conductance values cause
uncertainty in the forward field calculation. Thus, we can takeadvantage of the probabilistic framework of the proposed modelby including a hyperparameter over the forward field F thatwould enable the algorithm to make small changes to thecomputed F, using the data to “learn” an improved estimate of theforward field.
One could further relax assumptions about the source andinterference factors. They could be modeled by non-Gaussian
distributions, although care should be taken in choice of distribution
so that they could still be solved for analytically. The temporalprogression of the factors could be explicitly modeled, such as using
an auto-regressive model or Kalman filter.
The scanning framework of solving for one voxel's activity at atime can be relaxed. Spatial smoothness or spatial priors from othermodalities, such as structural or functional MRI, could beincorporated. Furthermore, one is not limited to sn in a single
voxel; the above formulation holds for any P arbitrarily chosen
dipole components, no matter which voxels they belong to, and forany value of P. Of course, as P increases the inferred value of Φ
becomes less accurate, and one might choose to restrict it to adiagonal or block-diagonal form.
Related methods for MEG and EEG source reconstruction
Source localization algorithms can be broadly classified asparametric or tomographic, making different assumptions to
overcome the ill-posed inverse problem. Parametric methods,
including equivalent current dipole (ECD) fitting techniques,assume knowledge about the number of sources and theirapproximate locations. A single dipolar source can be localized
well, but ECD techniques poorly describe multiple sources orsources with large spatial extent. Parametric methods could also
involve multipole, line source and/or patch reconstructions(Schmidt et al., 1999; Yetik et al., 2005; Jerbi et al., 2002). Many
parametric methods involve computation of noise and or signal
covariances from data where several Bayesian approaches havebeen reported (Mosher and Leahy, 1998; Jun et al., 2005, 2006;
Huizenga et al., 2002; Waldorp et al., 2002; de Munck et al., 2004;
Bijma et al., 2005).
Tomographic methods reconstruct an estimate of source activity
at grid points across the whole brain, either simultaneously or
through scanning. Tomographic methods could also be classified
as adaptive or non-adaptive. Of many tomographic algorithms, theadaptive beamformer has been shown to have the best spatial
resolution and zero localization bias (Sekihara et al., 2002, 2005;
Küçükaltun-Yildirim et al., 2006); however, it is sensitive totemporally correlated sources. Several tomographic algorithms use
a Bayesian framework to explicitly incorporate prior assumptions
about source configuration in the model (Sato et al., 2004; Phillips

et al., 2005; Mattout et al., 2006). Wipf et al. (2007) provide a
unifying framework relating methods through their choice in
specification of a neural source prior, whether fixed or learned
from the data.
In contrast to many of these methods which reconstruct wholebrain source activity as MAP estimates, SAKETINI is a scanning
based method that reconstructs a full posterior distribution of thesource activity at each scanning voxel. The closest algorithm
related to this approach is Dogandzic and Nehorai (2000),
however, their algorithm does not make use of latent variable
modeling and assumes that the temporal basis functions of thesource activity at a scan location are known. SAKETINI does notmake this assumption. Furthermore, instead of estimating a fullrank noise covariance, SAKETINI uses variational Bayesian factor
analysis to estimate a dimension-reducing mixing matrix withnumber of columns much less than number of rows. Thus, morerobust estimates of noise covariances can be made with fewer
parameters to estimate.
Jun et al. (2005) create a multi-dipole model that explicitly
computes the posterior for number of dipoles, location andorientation. They use prior terms for the temporal correlations of
the sources, chosen based on the data. A prior on location ororientation can include information from fMRI or structural MRI,
although they just use a uniform prior. The prior on number ofdipoles roughly scales with the inverse of the number. They
acknowledge the problem that the noise covariance has many
parameters to estimate, so they use a Wishart prior distribution andestimate the noise posterior. They also integrate over the timecourses to obtain the posterior for location, orientation and numberof dipoles.
In this paper, we have not included any specific spatial priorinformation, although it certainly can be incorporated. SAKETINI
only estimates one dipole at a time, by scanning through the voxel
grid, thus estimation of the number of dipoles is not explicitly
performed. Thresholding of the likelihood map can be viewed as aposterior probability map thresholding procedure.
Conclusion
We have described a novel probabilistic algorithm which performs source localization while robust to interference anddemonstrated its superior performance over standard methods ina variety of simulations and real datasets. The model takes
advantage of knowledge of when sources of interest are notoccurring (such as in the pre-stimulus period of a evoked response
paradigm). It learns the statistical structure of the interference
sources from the pre-stimulus period and then can suppress thesesignals in the post-stimulus period. Set in a probabilistic framework, the use of prior probabilities and hyperparameters assists inthe computation of unknown quantities. The flexibility of themodel allows for many possible extensions.
Acknowledgments
We would like to thank Kenneth Hild and Ben Inglis for theirhelpful discussions and comments on the manuscript, Sarang Dalal
for helping with NUTMEG programming, Mary Mantle, Anne
Findlay, and Susanne Honma for helping with data collection, andthe anonymous reviewers for their helpful comments on the
manuscript. This work was supported by NIH grants R01
NS44590, DC4855 and DC6435.

J.M. Zumer et al. / NeuroImage 37 (2007) 102–115

References
Attias, H., 1999. Inferring parameters and structure of latent variable
models by variational Bayes. Proc. 15th Conf. Uncertainty Artif. Intell.,pp. 21–30.
Beal, M.J., May 2003. Variational algorithms for approximate Bayesianinference. PhD thesis, University of London.
Beal, M.J., Jojic, N., Attias, H., 2003. A graphical model for audiovixual
object tracking. IEEE Trans. Pattern Anal. Mach. Intell. 25, 828–836.
Bijma, F., de Munck, J.C., Heethaar, R.M., 2005. The spatiotemporal MEG
covariance matrix modeled as a sum of Kronecker products. NeuroImage 27, 402–415.
Brookes, M.J., Stevenson, C.M., Barnes, G.R., Hillebrand, A., Simpson,
M.I.G., Francis, S.T., Morris, P.G., 2007. Beamformer reconstructionof correlated sources using a modified source model. NeuroImage 34,
1454–1465.
Bunch, P., Hamilton, J., Sanderson, G., Simmons, A., 1978. A free response
approach to the measurement and characterization of radiographicobserver performance. J. Appl. Photogr. Eng. 4, 166–172.
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., Sekihara, K., Nagarajan, S.S., 2006. Modified beamformers forcoherent source region suppression. IEEE Trans. Biomed. Eng. 53,
1357–1363.
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.
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,
2123–2128.
Dogandzic, A., Nehorai, A., 2000. Estimating evoked dipole responses in
unknown spatially correlated noise with EEG/MEG arrays. IEEE Trans.
Signal 13–25.
Gelman, A., Rubin, D.B., 1996. Markov chain Monte Carlo methods inbiostatistics. Stat. Methods Med. Res. 5, 339–355.
Ghahramani, Z., Beal, M., 2001. Graphical models and variational methods.
In: Opper, M., Saad, D. (Eds.), Advanced Mean Field Methods—Theory
and Practice. MIT Press.
Hamalainen, 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.
Huizenga, H.M., de Munck, J.C., Waldorp, L.J., Grasman, R.P.P.P., 2002.
Spatiotemporal EEG/MEG source analysis based on a parametric noise
covariance model. IEEE Trans. Biomed. Eng. 49, 533–539.
Ikeda, S., Toyama, K., 2000. Independent component analysis for noisy
data–MEG data analysis. Neural Netw. 13, 1063–1074.
Jerbi, K., Mosher, J.C., Baillet, S., Leahy, R.M., 2002. On MEG forward
modelling using multipolar expansions. Phys. Med. Biol. 47, 523–555.
Jordan, M.I. (Ed.), 1999. Learning in Graphical Models. MIT Press,
Cambridge, MA, USA.
Jordan, M., Ghaharamani, Z., Jaakkola, T., Saul, L., 1999. An introduction
to variational methods for graphical models. In: Jordan, M. (Ed.),
Learning in Graphical Models. MIT Press, Cambridge, MA, USA.
Jun, S.C., George, J.S., Paré-Blagoev, J., Plis, S.M., Ranken, D.M., Schmidt,
D.M., Wood, C.C., 2005. Spatiotemporal Bayesian inference dipole
analysis for MEG neuroimaging data. NeuroImage 28, 84–98.
Jun, S.C., George, J.S., Plis, S.M., Ranken, D.M., Schmidt, D.M., Wood,
C.C., 2006. Improving source detection and separation in a spatiotemporal Bayesian inference dipole analysis. Phys. Med. Biol. 51,
2395–2414.
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, 163–178.

115

Küçükaltun-Yildirim, E., Pantazis, D., Leahy, R.M., 2006. Task-based
comparison of inverse methods in magnetoencephalography. IEEE
Trans. Biomed. Eng. 53, 1783–1793.
Limpiti, T., Van Veen, B.D., Wakai, R.T., 2006. Cortical patch basis model
for spatially extended neural activity. IEEE Trans. Biomed. Eng. 53,
1740–1754.
Mattout, J., Phillips, C., Penny, W.D., Rugg, M.D., Friston, K.J., 2006. MEG
source localization under multiple constraints: an extended Bayesianframework. NeuroImage 30, 753–767.
Mosher, J.C., Leahy, R.M., 1998. Recursive MUSIC: a framework for EEG
and MEG source localization. IEEE Trans. Biomed. Eng. 45,
1342–1354.
Nagarajan, S., Attias, H., Hild, K., Sekihara, K., 2005. Stimulus evokedindependent factor analysis of MEG data with large background activity.
Adv. Neural Inf. Process. 18.
Nagarajan, S.S., Attias, H.T., Hild, K.E., Sekihara, K., 2006. A graphicalmodel for estimating stimulus-evoked brain responses from magnetoencephalography data with large background brain activity. NeuroImage
30, 400–416.
Neal, R.M. (Ed.), 1996. Bayesian Learning for Neural Networks. Springer,
New York.
Pascual-Marqui, R.D., 2002. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Meth. Find. Exp.
Clin. Pharmacol. 24 (Suppl. D), 5–12.
Phillips, C., Mattout, J., Rugg, M.D., Maquet, P., Friston, K.J., 2005. An
empirical Bayesian solution to the source reconstruction problem in
EEG. NeuroImage 24, 997–991011.
Sarvas, J., 1987. Basic mathematical and electromagnetic concepts of thebiomagnetic inverse problem. Phys. Med. Biol. 32, 11–22.
Sato, M.-A., Yoshioka, T., Kajihara, S., Toyama, K., Goda, N., Doya, K.,
Kawato, M., 2004. Hierarchical Bayesian estimation for MEG inverse
problem. NeuroImage 23, 806–826.
Schmidt, D.M., George, J.S., Wood, C.C., 1999. Bayesian inference
applied to the electromagnetic inverse problem. Hum. Brain Mapp. 7,
195–212.
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., Poeppel, D., Marantz, A., 2002. Performance ofan 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., Sahani, M., Nagarajan, S., 2005. Localization bias and spatial
resolution of adaptive and non-adaptive spatial filters for MEG sourcereconstruction. NeuroImage 25, 1056–1067.
Snodgrass, J.G., Corwin, J., 1988. Pragmatics of measuring recognition
memory: applications to dementia and amnesia. J. Exp. Psychol. Gen.
117, 34–50.
Vrba, J., Robinson, S.E., McCubbin, J., 2004. How many channels areneeded for MEG? Neurol. Clin. Neurophysiol. 2004, 99.
Waldorp, L.J., Huizenga, H.M., Grasman, R.P.P.P., Böcker, K.B.E., de
Munck, J.C., Molenaar, P.C.M., 2002. Model selection in electromagnetic source analysis with an application to VEFs. IEEE Trans. Biomed.
Eng. 49, 1121–1129.
Wipf, D., Ramirez, R., Palmer, J., Makeig, S., Rao, B., 2007. Analysis of
empirical Bayesian methods for neuroelectromagnetic source localization. In: Schölkopf, B., Platt, J., Hoffman, T. (Eds.), Adv. Neural Inf.
Process. Syst., vol. 19. MIT Press, Cambridge, MA.
Yetik, I.S., Nehorai, A., Muravchik, C.H., Haueisen, J., 2005. Line-source
modeling and estimation with magnetoencephalography. IEEE Trans.
Biomed. Eng. 52, 839–851.
Zumer, J., Attias, H., Sekihara, K., Nagarajan, S., 2007. A probabilistic
algorithm integrating source localization and noise suppression for MEG
and EEG data. In: Schölkopf, B., Platt, J., Hoffman, T. (Eds.), Adv.
Neural Inf. Process. Syst., vol. 19. MIT Press, Cambridge, MA.