www.elsevier.com/locate/ynimg
NeuroImage 41 (2008) 924 – 940

Probabilistic algorithms for MEG/EEG source reconstruction using
temporal basis functions learned from data ☆
Johanna M. Zumer,a,b Hagai T. Attias,c Kensuke Sekihara,d and Srikantan S. Nagarajana,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., P.O. Box 475608, San Francisco, CA 94147, USAd
Department of Systems Design and Engineering, Tokyo Metropolitan University, Tokyo, 191-0065 Japanb

Received 18 November 2007; revised 5 February 2008; accepted 11 February 2008
Available online 20 March 2008

We present two related probabilistic methods for neural source reconstruction from MEG/EEG data that reduce effects of interference, noise,and correlated sources. Both methods localize source activity using alinear mixture of temporal basis functions (TBFs) learned from the data.
In contrast to existing methods that use predetermined TBFs, we compute TBFs from data using a graphical factor analysis based model[Nagarajan, S.S., Attias, H.T., Hild, K.E., Sekihara, K., 2007a. A
probabilistic algorithm for robust interference suppression in bioelectromagnetic sensor data. Stat Med 26, 3886–3910], which separates
evoked or event-related source activity from ongoing spontaneousbackground brain activity. Both algorithms compute an optimal weighting of these TBFs at each voxel to provide a spatiotemporal map ofactivity across the brain and a source image map from the likelihood of adipole source at each voxel. We explicitly model, with two different
robust parameterizations, the contribution from signals outside a voxelof interest. The two models differ in a trade-off of computational speed
versus accuracy of learning the unknown interference contributions.
Performance in simulations and real data, both with large noise andinterference and/or correlated sources, demonstrates significant improvement over existing source localization methods.
© 2008 Elsevier Inc. All rights reserved.
Keywords: Biomagnetism; Magnetoencephalography (MEG); Electroencephalography (EEG); Inverse problems; Bayesian inference; Denoising☆

This work was supported by NIH grants R01 NS44590, DC4855 andDC6435.
⁎ 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 © 2008 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2008.02.006

Introduction
Magnetoencephalography (MEG) and electroencephalography(EEG) are popular methods for providing the spatiotemporal characteristics of human neural activity to both researchers and clinicians. Both techniques record the effects of neural activity at thescalp with millisecond precision. The increasing availability ofwhole-head MEG/EEG sensor arrays allows for higher-resolution
spatiotemporal reconstruction of neural activity, thus increasing thedemand for improved methods for source reconstruction.
Many sources of noise interfere with true signals in the MEG/
EEG data, affecting all existing inverse method algorithms. Thermal
or electrical noise is present at the MEG or EEG sensors themselves.
Background room interference such as from powerlines and electronic equipment can be problematic. Biological noise such as heartbeat, eyeblink, or other muscle artifact can also be present. Ongoingbrain activity itself, including the drowsy-state alpha (~10 Hz)
rhythm, can drown out evoked brain sources. Finally, many localization algorithms have difficulty in separating neural sources ofinterest that have temporally overlapping activity.
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 answered, and is prohibitive for examining processes such as learningthat can occur over just a few trials. Obtaining sufficient trials for
successfully averaging out noise is time-consuming and therefore
difficult for a subject or patient to hold still or pay attention throughthe duration of the experiment.
Algorithms proposed in this paper use a probabilistic graphicalmodel framework: a general tool for learning unknown, underlying
variables from observed sensor data. The graphical model depicts
probabilistic dependencies between nodes, which include the observed data, computed lead field, unobserved evoked and interference factors, and sensor noise. Many inference algorithms exist toJ.M. Zumer et al. / NeuroImage 41 (2008) 924–940

estimate the unknown quantities given the data and model. We haverecently shown that this approach is effective for interference suppression, source separation and source localization of MEG data
(Nagarajan et al., 2007a; Zumer et al., 2007).
The source reconstruction framework proposed in this paper isreferred to as Neurodynamic Stimulus-Evoked Factor Analysis Localization (NSEFALoc), which first uses a separate graphical model
called Stimulus-Evoked Factor Analysis (SEFA) to estimate temporal basis functions and then finds the best linear mixture (spatial
weighting) of these basis functions at each source voxel. NSEFALoc1 models the activity outside a particular voxel by a full-rank covariance matrix and estimates unknown quantities by maximizing
the likelihood. NSEFALoc2 parameterizes activity outside the voxelof interest as a linear mixture of a set of unknown Gaussian factors
plus Gaussian sensor noise and estimates all unknown quantities
using a Variational Bayesian Expectation-Maximization (VB-EM)
algorithm (Attias, 1999; Ghahramani and Beal, 2001). Both techniques create an image of brain activity by scanning the brain,inferring the models from sensor data, and using them to computethe maximized likelihood of the data with the best set of parametersat each voxel, creating a spatial map to indicate the most likely
locations of sources.
It is clear that improved performance for noisy data with correlated sources is a desirable trait for a new source reconstruction

925

method, especially since some methods such as minimum varianceadaptive beamforming (MVAB) is known to have reduced performance when at least two sources are highly correlated (Sekihara
et al., 2002a). The simulations and real data tested here illustratethese issues and demonstrate improved performance of NSEFALocover existing methods. In simulations, several parameters such as
location of sources, rotating or fixed dipoles, SNR, and type ofbackground noise were varied. The effect of number of sensors andtimepoints (total data available) was also tested. Finally, robustness
to choice of number of basis functions or factors by the user isshown. Furthermore, performance of all methods is compared using
some real-data examples from an auditory-evoked MEG data set anda low-SNR somatosensory MEG data set.
An initial report of this method was presented in Zumer et al.
(2006). This current paper expands on the mathematical details andprovides a more thorough analysis of performance in both simulations and real data, in comparison to established methods of MVAB
(Sekihara et al., 2002b) and sLORETA (Pascual-Marqui, 2002).
TheoryIn both NSEFALoc models, we assume the source activity is alinear combination of J × N temporal basis functions Φ computedfrom the data, spatially weighted at each voxel r by a Q × J dipole

Fig. 1. Graphical models for NSEFALoc1 and NSEFALoc2. Noisy sensor data are first processed by SEFA to determine the denoised temporal basis functions Φ,of which a linear mixture can produce any localized evoked source. These TBFs are then input as fixed bases to both NSEFALoc1 and NSEFALoc2, which
estimate the spatial weighting G of these TBFs for each voxel. The likelihood map can be displayed, and the source estimate at its spatial peaks can be plotted. In
each graphical model, quantities inside the large square are variables dependent on time while quantities outside are parameters/hyperparameters independent oftime. Directed arrows between nodes indicate a probabilistic dependence. Square nodes are known (observed or computed) while circles nodes are unknown. The
relative amounts of dashes/dots for each circle or square indicate groupings of nodes.

926

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

mixing matrix Gr. We compute the maximum likelihood at eachvoxel; the spatial peaks of this likelihood map correspond to themost likely source locations.
Fig. 1 depicts the graphical models for the processing steps ofboth NSEFALoc models. SEFA (top middle) is a separate modelthat is first run as a preliminary step on the data Y prior to either
NSEFALoc algorithm in order to learn the denoised evoked factorsΦ (top right) to be used at temporal basis functions. In SEFA,
evoked brain activity, biological noise, other room interference,and sensor noise all contribute to the measured sensor data (topleft). The second step is to run either or both NSEFALoc models,
shown in the second row of Fig. 1. Both NSEFALoc models use
the averaged sensor data and the temporal basis functions Φ as
known/fixed quantities. Finally, both models output a likelihood
map indicating the location of sources as well as the source timecourse estimates.
The mathematical notation used throughout this paper is asfollows. Matrices are in bold upper case, vectors are in bold lower
case (e.g. the nth column of the matrix Φ is / n, or a vector such as
the hyperparameter α), and scalars are in non-bold lower case (e.g.
the element from the kth row and lth column of matrix A is akl).
Non-bold upper case Roman letters are used to denote the dimension of matrices or vectors, such as K number of MEG sensors.

et al. (2007a) describe the one-stage model. Here, SEFA is computedby using a two-stage procedure, where B, U, and Λ S are first learned
from just the pre-stimulus data alone. Then, B and ΛS are held fixedand C, Φ, and U are computed using just the post-stimulus data. Notethat U needs to be recomputed for the post-stimulus period, since theprojection from the data to the noisy source space is defined by B,
which is fixed, but the actual realization of noise strength in the poststimulus changes from time point to time point.
We define  P  P P P fn
kn ¼; Ω¼ C B ; Ω¼ C B ;ð2Þun

Computing temporal basis functions from data using SEFANSEFALoc1

We assume that the neural activity at all possible source locations
can be described as a linear combination of temporal basis functions,which we estimate using SEFA (Nagarajan et al., 2007a). SEFA uses
the computational framework of Variational Bayesian Factor analysis (VBFA), but includes the additional concept of how MEG/EEG
data are collected in order to further separate out noise components.
In the stimulus-evoked paradigm, some baseline control data arecollected for several hundred milliseconds during the pre-stimulus
period, then a stimulus occurs, evoking a neural response in the poststimulus period.
The key idea of SEFA is that background activity, such as
ongoing brain activity unrelated to the stimulus, other biologicalnoise such as eye blinks and heartbeat, other room noise, andsensor noise, will be present in both pre-stimulus and post-stimulus
periods. However, only the evoked brain sources of interest will bepresent in post-stimulus period alone and not in the pre-stimulus
period. Note that this assumption is valid for the evoked response
paradigm, but not the event-related synchronization/desynchronization analysis.
The data Y is partitioned into pre- and post-stimulus sections as:

The NSEFALoc1 model and its solution are related to that proposed by Dogandzic and Nehorai (2000) and Baryshnikov et al.
(2004). NSEFALoc1 differs from their work by precomputing the–basis functions /n from the data using the estimated / n from SEFA.
These SEFA estimates are preferred since interference and noise
sources have been removed, and the spectral content and statisticalproperties have not been restricted. NSEFALoc1 also differs from theabove methods by placing a Wishart prior distribution on the fullrank precision matrix, which assists learning many unknown quantities from potentially few data points. The NSEFALoc1 generativemodel for the K × 1 sensor data yn is:yn ¼ Bun þ vn
yn ¼ Cfn þ Bun þ vn

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

ð1Þ

Time ranges from −Npre:0:Npost − 1 where Npre (Npost) indicatesthe number of time samples in the pre- (post-) stimulus period. The
K × M matrix B and the M × 1 vector u n represent the backgroundmixing matrix and background factors, respectively. The K × L
matrix C and L × 1 vector / n are the evoked mixing matrix andevoked factors (temporal basis functions), respectively. The sensor
noise term v n is described by diagonal precision matrix ΛS, where
the subscript S indicates the Λ learned from SEFA. The quantities
Φ, U, and V are the matrices for all time points for / n, u n, and v n.
The details of the model are described in Appendix A. The update
rules are listed explicitly here again, since those listed in Nagarajan

The main update equation for / n from the post-stimulus data
is:


P
qðkn jyn Þ¼ N kn jP
kn ;Γ ;Pkn ¼Γ1 Ω T ΛS yn ;PPP

Γ ¼ ΩT ΛS Ω þ I ¼ Ω T ΛS ΩþKΨ1 þ I 10P P PTΨC 0KþI
B C ΛS C B þ
0
0
¼ @PT
B

Þ

yn ¼ Fr Gr fn þ wrn

ð3Þ
ð4Þ

ð5Þ

Both NSEFALoc models are based on a physical description ofneural activity, in which brain sources are modeled by current
dipoles. For a given volume conductor forward model, the K × Q
forward lead field matrix Fr represents the physical relationshipbetween a dipole at voxel r and its influence on sensor k = 1:K(Sarvas, 1987). In the most general case, including EEG data, Q = 3
for all three possible directions of coordinate bases of a source
dipole. In the case of the single-shell sphere as commonly used in
MEG, the radial component of source dipoles contributes nothing
to MEG sensors, thus Q = 2. If there is knowledge of subjectspecific cortical anatomy, the source may be constrained to beperpendicular to the gray matter surface, thus Q = 1. Throughout
the rest of this paper, the single-shell model with Q = 2 is used forboth simulations and real data from MEG, although these methods
could be easily extended to a multisphere model for MEG with
Q = 3 or to EEG data with an appropriate forward model taking
tissue conductivities into account.
The noise wrn is modeled by zero-mean Gaussian distribution
with a K × K precision matrix Λ1 which is full-rank (not diagonal
like in SEFA above), and the subscript 1 indicates the Λ learned in
the NSEFALoc1 model. Both the parameters G and Λ1 are unknown.

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

927

For a large number of sensors K, the precision matrix becomes quite
large and difficult to infer accurately from the data. It may also
become ill-conditioned. Hence, a prior probability using a Wishart
distribution is used for Λ1:

The following conditional probabilities complete specificationof the model:pðΛ1 Þ ¼ W ðΛ1 jn; Σ0 Þ~jΛ1 jm=2 e2TrðΣ0 Λ1 Þpðxn Þ ¼ N ðxn j0; IÞ; pðvn Þ ¼ N ðvn j0; Λ2 Þ pðAÞ¼ p akj ¼ N akj j0; ðk2 Þk αj1

ð6Þwhere Σ0 and ν are hyperparameters. A Wishart distribution isrelated to a multivariate Γ distribution.
The estimates Λ̂1 and Ĝ are the value of each that maximizes
the likelihood L specific to each scanned voxel, given in Eq. (B1).
Their derivation is described in Appendix B and the final results
are given here. Initially solving for Λ̂−1 1 gives:
1

Λ̂1 ¼

1
ðRYY  FGRΦY þ Σ0 Þ
N þm

ð7Þ

Solving for G r gives
ð8Þ

N þm
logjΛr1 j þ const:2

ð9Þ

whose spatial peaks correspond to the most likely source locations.
The data and factor covariance matrices referred to above are:RYY ¼

yn yTn ;

n¼1

RYΦ ¼

N
X

yn fTn ;

RΦΦ ¼

n¼1

N
X

fn fTn

ð10Þ

n¼1

ð11Þ

where ΓΦ− 1 indicates the first set of rows only corresponding onlyto Φ, but all columns, and Λ S is from Eq. (A8).
NSEFALoc2–
NSEFALoc2 also uses the TBFs Φ estimated from SEFA described above. In contrast to NSEFALoc1, the contributions to poststimulus sensor measurements not arising from a dipole source at thevoxel r are now more explicitly modeled in NSEFALoc2. The J × 1
unknown interference factors x\rn correspond to activity in all voxels
excluding r and A\r is a K × J unknown mixing matrix, where \rmeans corresponding to activity not at r. The sensor noise has
unknown diagonal precision Λ2, where the subscript 2 indicates theΛ learned in the NSEFALoc2 model. The corresponding generativemodel for the sensor data Y is:ryn ¼ Fr Gr fn þ A⧹r x⧹rn þ vn



kj

kj

ð12Þ

ð14Þ

Notice that in place of the (K + K) / 2 elements of the full-rank
precision matrix Λ1 in NSEFALoc1, now just the KJ + K elements of
A and diagonal Λ2 need to be inferred from the data. Since typicallyJ≪K (J b ~10 and at UCSF K = 275), NSEFALoc2 has significantly
less parameters and can thus be inferred more accurately.
We again use a VB-EM algorithm to infer the unknown quantities from the data and the derivation is given in Appendix C. The
VBE-step updates for the variables are:
 P pðxn jyn Þ ¼ N xn j xn ;Γ
P

P

ð15Þ

Γ ¼A Λ2 A þ KΨ þ I
TIn the VBM-step, the full posterior over A is found by finding
the q(A|Y) that best approximates p(A|Y) and the MAP estimatesof the parameters G and Λ2 and hyperparameter α are found. The
posterior distribution of A is thus: Pqðak jYÞ; qðak jYÞ ¼ N ak j ak ;ðk2 Þk ΨqðAjYÞ ¼kP
ð16Þ
A ¼ðRYX  FGRΦX ÞΨ1
Ψ ¼ RXX þ αj

The MAP estimate of G is



1
G ¼ FT Λ2 F FT Λ2 RYΦ  RYX Ψ−1 RXΦ1 RΦΦ −RΦX Ψ1 RXΦP

The source estimate from both NSEFALoc1 and NSEFALoc2 isgiven by G rΦ. For NSEFALoc1, using Eqs. (8), (10), and (3), thesource estimate per voxel r is1
1 P Tŝ rn ¼ FrT S1 Fr FrT S1 RYΦ R1ΦΦ ΓΦ Ω ΛS yn



jp a ;

xn ¼Γ1 ATΛ2 ðyn  F Gfn Þ

Since the expression for G is now known, this value can beplugged into Eq. (7) to find Λ1. The maximized likelihood is then:

N
X

ð13Þ

2

P


1
G ¼ FT S1 F FT S1 RYΦ R1ΦΦ

1 
RYY  RYΦ R1
S¼
RΦΦ ΦY þ Σ0
N þm
r

Lr ¼

pðyn jxn ; A; Λ2 Þ ¼ N ðyn jFGfn þ Axn ; Λ2 Þ

ð17Þ

Solving for Λ2 and α in NSEFALoc2 is very similar to solvingfor Λ S and χ in SEFA, by letting y n′ = y n−FG/n. Then, take thederivative of the free energy F w.r.t. Λ2 (or α) to obtain a similarsolution:
h

i1
P
Λ2 ¼ N diag RY VY V  ARXY Vð18Þ1 PT Pα1 ¼ diagA Λ2 AþΨAK
The expressions RYX, RY′X , RΦX, and RXX represent theposterior covariance between the two subscripts, similar to matrices previously defined. The maximized likelihood function for
NSEFALoc2 is the following, where the dependency on voxel
location is made explicit:
Lr ¼

NjΛr =2pj K
log 2 r þ logjαr Ψr j
2
jΓ j
2
N

T
Pr
Pr
1X
P
P
yn  FrG fn Λr2 yn FrG fn  x Tn Γr xn

2 n¼1

ð

ð

Þ

Þ
ð19Þ

The source estimate for NSEFALoc2, using Eqs. (17) and (A9) is:



1
ŝrn ¼ FT Λ2 F FT Λ2 RYΦ −RYX Ψ1 RXΦ


PT
1
 RΦΦ  RΦX Ψ1 RXΦ Γ1Φ Ω ΛS yn

ð20Þ

928

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

Methods: simulations, performance metrics, and real data
Simulation setup
The construction of simulated data sets and performance metrics
was similar or identical to that described in Zumer et al. (2007).
Simulations were created using a variety of realistic source configurations reconstructed on a 5 mm voxel grid. A single-shell spherical volume conductor model for MEG data was used to calculate
the forward lead field (Sarvas, 1987). While EEG models and data
were not tested here, use of an appropriate forward model for EEG
would make NSEFALoc amenable to EEG data as well.
Simulations and real data were analyzed 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-ofinterest, computation of forward field, filtering, and other denoising
preprocessing methods, as well as a variety of source reconstruction
methods, including MVAB (Sekihara et al., 2002b), sLORETA(Pascual-Marqui, 2002), SAKETINI (Zumer et al., 2007), timefrequency methods (Dalal et al., 2007), and now NSEFALoc.
Gaussian-damped sinusoidal time courses at specific locations
inside a voxel grid are based on realistic head geometry. Sourceswere set to be active only during a post-stimulus period, which
always composed 62.5% of the total data available, while theremaining 37.5% were pre-stimulus data. Typically 700 total data
points were used, unless specified otherwise.
In some simulations, only Gaussian sensor noise was added to theprojected simulated sources, termed the sensor noise only case. While
this type of simulation is common in simulation testing, this clearly
does not reflect true data which has interference sources somewhere
in source space contributing to covariance across sensors.
In another set of simulations, termed simulated interference cases,
background activity in source space was drawn from the Gaussian
distributions assumed by the model to simulate ongoing brain activity.
These background sources were placed in 30 random locations
throughout the brain voxel grid, active in both pre- and post-stimulus
periods. Their activity was projected onto the sensors and added to bothGaussian sensor noise and source activity. These simulated backgroundbrain 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, a final set of simulations was created termed real brain noise. Real MEG sensor data
were collected from a CTF MEG System with 275 axial gradiometers
while a human subject was alert but not performing tasks or receiving
stimuli. These background data thus include real sensor noise plus
real ongoing brain activity that could interfere with evoked sources
and add spatial correlation to the sensor data. Since throughout thiswork averaged data are used, these real data were binned into 100
trials of 700 data points each and averaged. The output signal to noise
ratio (SNR) and the corresponding output signal to noise-plus-interference ratio (SNIR) were varied. Output SNIR is calculated from theratio of the sensor data resulting from sources only to the sensor datafrom noise plus interference, as shown in the first equation of theResults section of Nagarajan et al. (2007b).
Localization accuracy of a single source
A single source was placed randomly within the voxel grid spaceand projected to the sensors, and all three types of noise (sensor noise

only, simulated interference, and real brain noise) were added to thesimulated sensor data, at four different levels of SNIR. Twenty
different realizations of random location were tested for each SNIR.
The localization error (in Euclidean distance) between the maximum
peak in the NSEFALoc1 and NSEFALoc2 likelihood maps, as wellas the MVAB and sLORETA power maps (sum of squares of poststimulus time points), and the true simulated source location wasmeasured.
Constructing simulations with three sources
Multiple active sources are more realistic than a single active
source. Moreover, any additional active source acts as interference
toward the ability to localize the first source. Several simulation
parameters were varied across different simulations. (i) Two different source configurations were used: one with three sources near
the surface as depicted in Fig. 3 and the other configuration withthree deeper sources. (ii) The orientation of the source was fixed inhalf the simulations and allowed to rotate over time in the otherhalf. (iii) The correlation of two of the three sources with eachother was set to be ρ = 0, ρ = 0.95, or ρ = 1; the third source wasalways uncorrelated with the other two sources. (iv) Each combination of parameters was tested for 10 different randomly generated
source time courses and source orientations. (v) In addition to thetrue source contribution to the sensor data, the three cases of sensor
noise only, simulated Gaussian interference or real brain noise were
tested. (vi) SNIR was set at 5 dB, 0 dB, or − 5 dB, with corresponding SNR for each case of 10 dB, 5 dB, and 0 dB. Thus, a totalof 1080 simulations were run using all combinations of simulation
parameters.
Model order selection and hyperparameters
The experimenter is faced with choosing the number of temporal basis functions for NSEFALoc methods (i.e. number offactors in SEFA) and the number of non-localized evoked factor X
in NSEFALoc2; the effect of this choice on performance should betested. The hyperparameters are expected to zero-out extra dimensions to a large extent, but this also should be tested. Finally, thechoice of SEFA for computing TBFs from the data should becompared to other options.
The choice of dimension for X (non-localized evoked factors)
in NSEFALoc2 was run with either 25 or 10 dimensions for A. The
inverse hyperparameter α− 1 over the mixing matrix A was normalized to the first hyperparameter for each of 1760 voxels examined. The localization result for both dimension choices was
examined.
The number of TBFs used in both NSEFALoc methods was
next tested and averaged over many simulations. In all simulations,
three sources are present with fixed dipole orientation. In half, twoof the three sources are perfectly correlated, while the other half ofthe simulations has uncorrelated sources; thus either two or three
TBFs are needed, respectively. Performance is characterized as
described in the next subsection and is compared to MVAB andsLORETA.
From Nagarajan et al. (2007a), SEFA seems to be a very good
way to obtain temporal basis functions of the denoised evoked
activity from real data when the true time course is not known.
However, their use as an input to the NSEFALoc class of models is
tested here. The simulations were the same as in Fig. 7 with threesources placed and three levels of correlation between two of theJ.M. Zumer et al. / NeuroImage 41 (2008) 924–940

three sources, and either rotating or fixed dipole orientation. The
use of SEFA to obtain TBFs was compared with using PCA andthe true time courses. The number of temporal basis functions washeld the same across the three different types of TBFs for eachsimulation (fewer were used when sources were known to becorrelated).
Performance evaluation
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, the ROC (receiveroperator characteristic) method was modified for brain imaging
results as suggested by Darvas et al. (2004), which is a measure of
hit rate versus false positive rate. The free-response ROC (FROC)
curve in particular allows for multiple hits per image (Bunch et al.,1978).
A local peak is defined here as a voxel that is greater in value thanits 26 three-dimensional neighbors. A hit is defined as a local peakthat is within a specified distance of the true location and above acertain threshold. A miss is defined as a true source location that has
no hit within the specified distance. A false positive is a local peakabove a certain threshold but farther than the specified distance froma true source location. A true negative is any voxel that is none of theabove.
FROC curves are generated by varying the threshold and allowable distance error, thus varying the trade-off of sensitivity andspecificity. The following distances were used as allowable localization error pofﬃﬃﬃa local peak
pﬃﬃﬃto a true locationpﬃﬃﬃ in order to be counted
as a hit: 54 3mm, 104 3mm, or 154 3mm. The threshold wasvaried to be 30%, 50%, 70%, or 90% of the maximum value in thewhole image. Thus, a hit rate (HR) and a false positive rate (FR)were recorded for each of the 12 combinations of threshold/error
for each of the 1080 simulations.
Since these HR versus FR points do not increase monotonically,as they would if threshold was the only criteria varied, we chose to
use the measure of A′ (similar to use in Zumer et al., 2007). A′ is a
way to approximate the area under the FROC curve for one HR/FR
point (Snodgrass and Corwin, 1988). The larger the area under the
FROC, the better the method is performing, since this means ahigher HR relative to FR for specified thresholds/localization errors.
For each simulation, the 12 computed A′ values were averaged togive one A′ value per simulation. The NSEFALoc1 and NSEFALoc2
likelihood maps were used as the spatial maps to test localization; thepower maps were used for MVAB and sLORETA.
Simulation: effects of number of sensors and time points
Previous studies have shown the advantage of sensor arrayswith larger number of channels (Hamalainen et al., 1993). Likewise, increased amount of data points across time usually leads toimproved estimation of unknown quantities. Therefore, the next setof simulations sought to determine how few sensors and how few
time points were needed to preserve performance.
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 the929

CTF system. Here, only a random subset of sensors was selected,using 150, 74, or 37 sensors. The numbers 74 and 37 were specifically chosen to correspond to the BTi commercial MEG system
previously installed in the UCSF lab until 2004.
To test the effect of the number of data points, the full set of 275
channels was used, but the available amount of data points was
reduced. All previous simulations have used 700 total data points,
where 62.5% were in the post-stimulus period. The ratio of data
points in the post-stimulus period was kept the same, but the total
number was reduced to 300, 200, 150, 100, or 50 time points.
Real MEG data
Several real data sets 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 gave written, informed consent toparticipate in each study, according to UCSF institutional review
board approval.
Auditory data sets were obtained by presenting 120 repetitions ofa 1 kHz tone binaurally to healthy subjects, 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, known to cause problems for the MVABTs
ability to localize the auditory sources.
We next examine a somatosensory data set 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 subjectTs 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 sourcereconstruction methods. We first applied NSEFALoc1, NSEFALoc2, MVAB, and sLORETA to the average of all 256 trials toassess the performance for the standard (high) SNR case. We thenapplied all three methods to the average of only the first 5 trials. To
further test if the performance was consistent across other sets ofjust 5-trial averages, we applied the three methods to the 5-trial
average of trials 6–10, 11–15, and 16–20. We then averaged theresults of these four different results. Any location found consistently will show up in the average.
Results
Single source localization
The mean localization errors for a simulated single source areshown in Fig. 2. Even at the lowest SNIR of − 5 dB, NSEFALoc1 and
NSEFALoc2 localized the source to within 5 mm error, which, forreal data, is on the order of the error due to coregistration of MEG
data with the subjectTs MRI. For all values of SNIR, NSEFALoc1
and NSEFALoc2 resulted in reduced error compared to MVAB andsLORETA. Errors for the sensor noise only case were not shownsince they were zero or essentially zero for all methods for all valuesof SNIR.
Examples of multiple sources, including correlated sources
Next, performance of the proposed models was tested for threesimultaneously active sources. Fig. 3 shows performance in two
examples each with either three uncorrelated sources (top half) or930

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

Fig. 2. Average localization error over 20 realizations of a randomly placed single dipole source. Background activity was either simulated interference or realdata. The standard error was typically 1 mm, not larger than 4 mm; error bars were omitted from the plot.

with 2 of 3 sources correlated (bottom half). All sources were fixedin orientation across time, with real brain noise added with SNIR of
5 dB. The only difference between the two examples in each half isthe random realization of the source time course. Two examples areshown to illustrate how a change just in temporal dynamics canaffect localization results. For each example, the log likelihood (orpower) map is above a grouping of plots showing the estimatedtime courses (gray) of the three sources overlaid onto the true timecourses (black). Note in the MVAB power maps that the threesources are labeled with a square, diamond, and circle and the timecourses plots are labeled accordingly. The bottom half of the plots
shows two square sources and a circle source, that indicating thetwo square sources are highly correlated.
The uncorrelated-source examples in Fig. 3 (top half) show that
all methods localize all three sources either perfectly or near
perfectly. The spatial peaks for MVAB are so focal they are hiddenby the square, diamond and circle symbols. The sLORETA powermap shows some difficulty in finding the lower left source perfectly, but does show a peak nearby. The NSEFALoc1 likelihood
map finds all three sources perfectly, though in Example 2, there isa possible false positive around (x = − 25, z = 5). The NSEFALoc2
likelihood map also finds all three sources perfectly, but also with apossible false positive around (x = − 30, z = 35). Both NSEFALoc1
and NSEFALoc2 display log-likelihood maps, which leads to theirincreased spatial spread relative to the MVAB power map, but ofcourse does not affect location of peaks.
The top half of Fig. 3 also shows all methodsT ability to estimate
the source time course. These examples show that sLORETA estimates the shape and amplitude very well for all sources. NSEFALoc2
also estimates the times courses well, although there is some crosstalk in Example 1, diamond source. NSEFALoc1 shows more severe
crosstalk errors. MVAB estimates the time courses reasonably well,although a slight misestimation of amplitude is seen.
The performance of all methods was further tested when two of
the three sources are highly correlated in time. The lower half ofFig. 3 shows the results from two simulation examples where twosquare sources are highly correlated in time (ρ = .95). The same real
brain noise was added at SNIR = 5 dB and all other aspects of thissimulation were the same as the uncorrelated case above, exceptthat the right square source time course was adjusted to correlatestrongly with the left square source time course. In all cases, theestimated time course plotted is the one extracted from the true
location, regardless of the localization map peak locations.
These examples illustrate the failure of MVAB for correlated
sources. The MVAB power map in Example 1 finds the uncorrelated
circle source, but largely mislocalizes the first square source and only
weakly finds the second square source. The reduction in power is seenin the first square time course plot. The MVAB power map for
Example 2 localizes all three sources within a reasonable error;
however, the amplitude of the peak location the two square correlated
sources is much reduced and might not be detected depending on thethreshold, which is also indicated by the large reduction in time series
amplitude. On the other hand, sLORETA is not in theory supposed to
be sensitive to correlated sources; sLORETA finds all three sources in
Example 1 (though one is weak and a center-of-the-head false positive
is of larger amplitude) but, in Example 2, fails to show distinct peaks
for the two correlated source locations. Despite these localization
issues, the time course estimation of sLORETA in both shape andamplitude was very accurate.
Overall, NSEFALoc2 localizes the sources and estimates source
time courses better than MVAB and sLORETA in these examples ofcorrelated sources in the bottom half of Fig. 3, while NSEFALoc1
performance is in-between. NSEFALoc1 localizes the sources wellin Example 1, but fails to find the first square source in Example 2,with a distant false positive instead. The time course estimatesby NSEFALoc1 suffer from quite a bit of crosstalk. Finally, the
NSEFALoc2 likelihood maps localize all three sources clearly with
only a slight localization error in the first square source. Furthermore,
the NSEFALoc2 time course estimation is quite good in both shape andamplitude, with much less crosstalk errors than NSEFALoc1 eventhough the same set of temporal basis functions was used.

Performance evaluation results
The performance of the proposed methods is now shown according to the metrics of A′ (area under ROC curve) and time course
estimation. Fig. 7(a) plots A′ for each method, for each value ofsource correlation and SNIR, and for all types of interference.
NSEFALoc1 and NSEFALoc2 both show A′ higher than MVAB andsLORETA. For the perfectly correlated source cases, NSEFALoc2
localizes sources best for all noise types.

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

931

Fig. 3. Performance of all methods in several example simulations. In the top two examples, all three true source locations, marked by square, diamond, andcircle, are uncorrelated with each other. On the bottom half, the true sources labeled with squares indicate location of true sources highly correlated with eachother, while the circle source is uncorrelated with the other two. While the source locations are the same for all examples, the time series are different for each.
Intensity of map corresponds to a normalized log-likelihood map for NSEFALoc1 and NSEFALoc2, and a normalized power map for MVAB and sLORETA.
Below the localization map for each example, black lines indicate simulated time series for each of the three source locations; gray lines indicate estimates of thesource time series at those three locations. The labels of squares, circles, or diamonds are included in each time series plot to indicate correspondence with thelocation on the map. The correlation of the true time course and the estimated time course is shown next to the symbol within each time series plot.

932

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

Fig. 4. Performance of NSEFALoc1 and NSEFALoc2 relative to MVAB and sLORETA for a variety of simulated data sets. Each data point is an average of 40
simulations, consisting of two different source locations and either a fixed or rotating source orientation. Standard errors were less than 0.05 for all points (notshown). (a) A measure of area under ROC curve A' is plotted in 9 subplots as a function of SNIR for sensor noise only, simulated and real brain interference(across columns), and for each of three source correlation values (across rows). See text for discussion of the A' metric. (b) The correlation of the estimated withthe true time course is plotted for each method.

The other main test of performance was the ability to estimate thesource time course. The estimated time courses for all methods wereobtained from the true source locations, regardless of whether theirrespective localization maps found that source as a hit. The correlationof the true time course with the estimated time course was computedfor each simulation and the averages are plotted in Fig. 7(b). In thesensor noise only case, sLORETA (dashed) estimates the time course
better than other methods regardless of source correlation, as previously understood for this method. However, when other sourcespace interference or real brain noise is added, this advantage ofsLORETA is lost. Instead, similar to the A′ results, NSEFALoc2
estimates the source time course the best when sources are perfectly
correlated in simulated interference and real brain noise cases.
Model order and basis function selection results
The ability of SEFA to learn the correct dimension of evoked
activity through the hyperparameters is demonstrated in Figs. 11 and
12 (Nagarajan et al., 2007b), and so is not examined further here.
Fig. 4 shows examples of NSEFALoc2 performance while thedimension (number) of non-localized evoked factors X is varied. The
main plots are of α− 1 (which controls number of evoked factors not atthe voxel of interest) normalized to the first hyperparameter (which isomitted from the plot). Each line within all plots is the value for each of1760 voxels analyzed. NSEFALoc2 was run with either 25 (left plot)or 10 (right plot) dimensions for A. From using 25 dimensions, it
seems that the inverse hyperparameters get close to zero after about 10
and that the extra dimensions are not contributing much. By using only
10, the values stay roughly the same but do change somewhat.

Furthermore, the localization results (inset in each plot) are roughly thesame yet some differences exist. The left inset (25 dimension result)shows all three sources with strong likelihood, but the lower source isblurred with the source above it. The right inset (10 dimension result)shows all three sources as distinct peaks although the lower source isweaker, and also has a possible false positive near it.
Fig. 5 compares NSEFALoc1 and NSEFALoc2 with three typesof temporal basis functions: the true source time sources, thoseobtained from SEFA (as the models were intended), or from PCA.
The performance metrics of A′ and time course estimation were used
to compare the choice of TBFs. In real brain noise and simulated
interference, using PCA to obtain TBFs resulted in the worst performance for both metrics of A′ and time course estimation. In
comparing SEFA with the true time courses, A′ is not affected, buttime course estimation is worse when using SEFA compared to thetrue; however, NSEFALoc2 with SEFA TBFs performs reasonably
close to the true TBFs (while NSEFALoc1 is considerably worse).
Finally, Fig. 6 shows the performance (through the A′ metric andtime course estimation) of all the methods as the number of dimensions
was varied, averaged over many simulations. For MVAB, the x-axis
represents the number of eigenvalues; for NSEFALoc1 and NSEFALoc2, it is the number of temporal basis functions, and it is meaninglessfor sLORETAwhose performance does not depend on such a parameter.
Overall, the lines for all methods are relatively flat, indicating
not too large of a dependence on the number of dimension reduction. The time course estimation in sensor-noise only shows theclear improvement in using at least three dimensions. The
correlated-source-case (bottom row) in interference or real brain
noise shows the clear advantage of NSEFALoc2 over MVAB fortime course estimation. Interestingly, the A′ metric gets worse forJ.M. Zumer et al. / NeuroImage 41 (2008) 924–940

933

Fig. 5. Plots of α- 1 hyperparameter for NSEFALoc2. Each line within all plots is the value for each of the 1760 voxels analyzed. The first hyperparameter isnormalized to one but is not shown. Inset in each plot is the localization result with the given number of dimension chosen, and symbols indicate correct location.
All three sources were rotating orientation and uncorrelated with each other, thus six independent time courses contributing to the sensor data. The left plot shows
results when the dimension of A was set to 10, while the right plot shows the dimension of A set to 25.

NSEFALoc2 in correlated sources in real brain noise as the number
of TBFs increases; this is probably due to incorrectly trying to fit
the extra components to the wrong location confused by thecorrelated sources.
Results as number of sensors and time points is varied
Fig. 8 shows simulation performance resulting from reduced
number of sensors. The top row is for SNIR = 0dB and the bottom
row is for SNIR = −10 dB. NSEFALoc1 and NSEFALoc2 did notshow any major degradation in performance for either A′ (leftcolumn) or time course estimation (right column) in the moderate
SNIR value of 0 dB. For the very noisy case of SNIR = − 10 dB, A′
begins to decline more with only 37 sensors; time course estimation
for the noisy SNIR = − 10 dB is poor for all number of sensors. In
contrast to the probabilistic methods, both MVAB and sLORETA
show decline in performance for both measures in the reduction from
275 to 150 sensors, but then plateaus for fewer sensors. In all cases of
150 sensors or fewer, both NSEFALoc methods outperform MVAB
and sLORETA for both metrics.

Fig. 6. Performance of NSEFALoc1 and NSEFALoc2 as a function of three types of temporal basis function used: the true source time sources, those obtained
from SEFA (as the models were intended), or from PCA. (a) A' metric for localization ability. (b) Time course estimation accuracy (similar to Fig. 7).

934

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

Fig. 7. Performance as a function of number of temporal basis functions for simulations with 3 fixed-orientation dipole sources. Source correlations of 0 and 1
were tested. Each data point is averaged over three values of SNIR (5, 0, - 5 dB). (a) A' and (b) correlation of estimated and true time course.

Fig. 9 shows the performance results of all methods with decreased
number of time points available. The top row is for SNIR = 0 dB
and the bottom row is for SNIR = −10 dB. The A′ results (left column)

show that both probabilistic methods outperform MVAB andsLORETA for all numbers of total data points; A′ performance
begins to decline for 150 or fewer data points. sLORETA is a non-Fig. 8. A′ and 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 41 (2008) 924–940

data-dependent method, thus the inverse weight is not affectedby number of time points available. The MVAB is dependent on thedata to provide an estimate of the data covariance matrix. Since
the simulations in both top and bottom rows are with relatively
high noise (SNIR = 0 dB and −10 dB, respectively), the data
covariance estimate might not change much with decreased data,
since it is already noisy (note the time course correlation does notreach above 0.5 for any number of data points tested at SNIR =−10 dB).
On the other hand, the time course estimation results (rightcolumn) show that NSEFALoc1, NSEFALoc2, and sLORETA donot show a decline in performance with fewer data points and that
all three methods generally perform equally well and better than
MVAB. This is most likely due to NSEFALoc methods not
requiring many data points in the first step of using SEFA to findtemporal basis functions; once the temporal bases have been found,less data are then needed for localization of these bases.
Somatosensory results
The left panel of Fig. 10(a) shows typical somatosensory-evoked
MEG data with the largest peak at 50 ms, expected to be coming fromprimary somatosensory cortex in the posterior wall of the central
sulcus. The next four panels of Fig. 10(a) show the localization
performance of NSEFALoc1, NSEFALoc2, MVAB and sLORETA.
All four methods accurately localize activity to the contralateral
primary somatosensory cortex. However, performance changes when
only 5 trials are used in the average. The left panel of Fig. 10(b) shows
the sensor data averaged over trials 1–5 of the same somatosensory
data set. The next four panels of Fig. 10(b) show errors in the935

localization in all methods. NSEFALoc1 and NSEFALoc2 showless error than MVAB and sLORETA, relative to the peak location
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, both NSEFALoc1 and NSEFALoc2
showed localization closest to the primary somatosensory cortex, as
shown in Fig. 10(c), whereas MVAB and sLORETA mislocalize
this source.
Auditory results
Fig. 11 shows localization results from all methods in four
different subjectsT AEF data sets. NSEFALoc1 finds activation inbilateral auditory cortex in 4/4 subjects, though extra peaks
appear in 3 subjects, and in Subject 3 the activation is too
superior. NSEFALoc2 finds activation in the bilateral auditory
cortex in 3/4 subjects, with extra peaks in only one subject, andan additional subject in which only the left auditory cortex isfound. The strongest peak in the MVAB power maps in all 4
subjects is (falsely) in the center of the head, while a weakeractivation on just the right side is seen in one subject. Finally,sLORETA finds bilateral auditory cortex in 2/4 subjects with extra
peaks in one of the two, and only one side of auditory cortex isfound in two other subjects. The sensor data for Subject 3 show
strong activation on the left side while very weak activation on theright, thus difficult to find for any method. While it is possible thatthe extra peaks seen in any of the methods are true sources coactivated with the primary auditory cortex, the sensor data do notgive a strong indication of extra sources, so most likely these extra
peaks are false positives. In general, NSEFALoc1 and NSEFALoc2

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

936

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

Fig. 10. 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 NSEFALoc1, NSEFALoc2, MVAB, and sLORETA. Row (a) shows performance of thethree 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 over other 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:
simulations, performance metrics, and real data for details. Crosshairs in localization maps show peak location within “active” voxels at the slice of the peaklocation, where the threshold for “active” was defined at 90% of the maximum for all maps.

found more correct source locations relative to less extraneous
peaks than MVAB and sLORETA.
Discussion
Two methods are introduced which localize stimulus-evoked
MEG/EEG sources and estimate their temporal activity in a probabilistic framework. Both model the sources as a linear combinationof denoised temporal basis functions derived from the data using avariational Bayesian factor analysis method. The methods have
reduced localization error relative to MVAB and sLORETA and are
not as hampered by correlated sources. Additionally, the number or
location of sources does not need to be specified, as in a standard
dipole fitting method. Thus, these methods have clear advantages
over current standard methods.
We showed results for MEG data only, although the equations canbe easily applied to EEG data with an appropriate lead field. Sources
can be constrained in location and orientation using the subjectTs
cortex defined by a structural MRI. Furthermore, NSEFALoc could
be modified to work with an extended lead field based on spatial
patch bases (Limpiti et al., 2006).
We have shown that the NSEFALoc models are not as sensitive
to temporally correlated sources as the standard formulation of
MVAB. However, it is possible to reduce the MVABTs 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).

As the number of MEG and EEG channels has increased in recent
years, the ability to accurately localize sources throughout the brain
has increased (Vrba et al., 2004). However, performing calculationsof high-dimensional data, such as inverting a data covariance matrix,
becomes more difficult and can lead to errors. Meanwhile, thedimensionality of the underlying neural activity remains the same.
Thus, many variations of PCA and ICA have been used on MEG/
EEG data for removal of noise/artifactual components as well as fordata dimension reduction (Jung et al., 2000; Ikeda and Toyama,
2000). Factor analysis also aims to reduce the dimensionality of thedata to a linear mixture of factors that best account for the data whileaccounting for noise at the sensor level. An extended version,
stimulus-evoked factor analysis, has been used here to partition thefactors that are event-related activity from the factors are background interference. All methods which perform dimension reduction need a criterion for choosing the reduction number. Using PCA,a plot of eigenvalues can often give a reasonable intuition for thedimension of “signal” in the data. ICA has no ordering of components. In the method proposed here, there are two variables affecting
model dimension: number of TBFs obtained from SEFA and thenon-localized evoked factors (X) in NSEFALoc2. While a user must
initially select a dimension for these terms, we showed that the use of
hyperparameters in the model provides robustness to this selection
by reducing the influence of unnecessary components.
For simple models with latent variables, the posterior distributionof a desired unknown variable can often be computed directly.
However, for more interesting and realistic models, the posterior is

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

often computationally intractable. In these cases, some approximation must be made. Since both SEFA and NSEFALoc2 models werecomputationally intractable as initially developed, we used a variational approximation for the joint posterior in both models. The main
alternative to variational methods is sampling methods, such as
Markov Chain Monte Carlo methods (Jun et al., 2005; Gelman and
Rubin, 1996), which extensively estimate points in the distribution.
MCMC is dependent on the sampled points and can be quite computationally costly. Nummenmaa et al. (2007a) do show advantagesof MCMC over variational methods when the posterior distribution is
not unimodal. However, the same researchers also show improvements of variational Bayesian methods over minimum norm methods
in real data (Nummenmaa et al., 2007b). Variational methods insteadchoose to factorize the joint distribution over factors and parametersassuming conditional independence of the factors and parameters,also termed the mean field approximation. Variational Bayesian
methods compute the posterior distribution that maximizes the freeenergy F, an approximation to the data likelihood L. This appro-

937

ximation is an equality when the approximate posterior q equals thetrue posterior p.
Several other uses of variational Bayesian methods for the MEG/
EEG inverse problem have been demonstrated. In general, they varyin how spatial priors, source covariance, and noise covariance aretreated, as well if they are a dipole or distributed source model. Satoet al. (2004) show how variational Bayesian inversion methods canbe used to improve MEG estimates with the inclusion of fMRI data.
Kiebel et al. (2008) use a variational Bayesian model for dipole
models; one benefit is to contrast competing models of number andtype of dipoles to overcome the usual problem of dipole models inchoosing the number of dipoles. Phillips et al. (2005) demonstrate a
distributed source model that uses multiple source priors and learns
their optimal weighting through hyperparameters. Friston et al.
(2008) extend this further by establishing a multiple prior formulation where any number of source prior covariances can be included, but are projected to sensor space, and their corresponding
hyperparameters prune which prior terms are relevant in sensor

Fig. 11. Performance of methods on real auditory-evoked MEG data sets from two healthy human subjects. NSEFALoc1 and NSEFALoc2 results are likelihood
maps and MVAB and sLORETA are power maps. The thresholds were set to portray each method optimally (i.e. including as many true sources as possible while
not including other areas).

938

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

space, thus avoiding large source-space matrices. Daunizeau andFriston (2007) use a variational inversion scheme to solve a multiscale model for MEG/EEG where the quantity and functional connectivity between mesostate sources are learned. Trujillo-Barreto et al.
(2008) have recently proposed a model that is similar to NSEFALoc inthat it includes a set of temporal basis functions to model the source
activity and accounts for the sensor noise and source noise separately
in a probabilistic graphical model; unknown quantities are also
learned through a VB-EM algorithm. Their method differs in several
ways from NSEFALoc. They demonstrate their method using
wavelet representation for TBFs; alternatively, SEFA could be used
to estimate the TBFs in their model. They estimate source activity at
all voxels at once rather than scanning each voxel at time.
NSEFALoc1 and NSEFALoc2 present a trade-off of computation
time and source estimation accuracy. Throughout the results presented
here, NSEFALoc2 tended to outperform NSEFALoc1. NSEFALoc1
estimates a full-rank noise covariance guided by a Wishart priordistribution in a single closed-form solution per voxel. NSEFALoc2, on
the other hand, learns more precisely the unknown interference sources
distant from the current voxel being scanned by learning an unknown
mixing matrix with a dimension smaller than the number of sensors.
Thus, more robust estimates of noise covariances can be made withfewer parameters to estimate, though convergence usually requires
about 20 EM iterations. These EM iterations require longer computation time: NSEFALoc1 computes estimates across a whole brain
volume of about 11,000 voxels in roughly 5 min, while NSEFALoc2
takes 110 min for the same reconstruction, roughly 0.6 s per voxel on astandard Linux personal computer with 2.0 GHz processor.
All methods which do not have a closed-form solution require theinitialization of the values to be iteratively updated. We have foundthat 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 that worked well in a few
test simulations, that set was used for all results shown. Since theclosed-form solution to NSEFALoc1 is easily obtained for eachscanned voxel, aspects of this result were used to initialize quantities
for NSEFALoc2, thus explaining some similarity in performance.
The NSEFALoc algorithms presented in this paper have somesimilarity to another algorithm, SAKETINI, recently proposed by us(Zumer et al., 2007). Both SAKETINI and NSEFALoc solve for
hidden evoked factors, use unknown mixing matrices to modelinterference sources, and take advantage of stimulus timing. However, SAKETINI does not use fixed temporal basis functions, butinstead learns hidden factors at each time point. Since the source
time course estimates from both NSEFALoc and SAKETINI are
effectively a weight matrix multiplying the sensor data, the temporalsmoothness of the source estimates is comparable to the sensor data,
possibly smoother due to noise removal. However, since the source
estimates from NSEFALoc are based on fixed temporal basis functions, additional smoothness could be imposed to these basis functions prior to source estimation; SAKETINI is not as amenable tothese modifications. The analysis of NSEFALoc is similar to thesimulations and real data that SAKETINI was tested on in Zumer et al.
(2007). A detailed comparison of performance between NSEFALoc
and SAKETINI is forthcoming and is beyond the scope of this paper.
In this work, no specific spatial prior information was used, although it certainly can be incorporated. NSEFALoc only estimates one
dipole at a time, by scanning through the voxel grid, thus estimation of
number of dipoles is not explicitly performed. The likelihood map canbe interpreted as a factorized map of the posterior probability of asource at each voxel. Thresholding of the likelihood map can be viewed

as a posterior probability map thresholding procedure. This posterior
probability map lends itself for statistical analyses across subjects andconditions, a topic that could be explored in future work.
Acknowledgments
The authors would like to thank Kenneth Hild and Ben Inglis
for helpful discussions including naming of the algorithm, Sarang
Dalal for help with NUTMEG programming, and Anne Findlay
and Susanne Honma for help with data collection.
Appendix A. Full set of update rules for SEFA estimates
Here, SEFA is computed by using a two-stage procedure to avoid
issues of identifiability between B and C, especially in cases of limited
pre-stimulus data. In the limit of no pre-stimulus data, B and C could
be concatenated as the same variable, as interference could not be
distinguished from evoked activity. However, even with sufficient
pre-stimulus data, if B and C are learned simultaneously, an evokedcomponent present only in the post-stimulus data could inadvertantly
be learned as a column of B due to the identifiability in this model.
To describe the full model in the Bayesian framework, priorprobability distributions are given to these quantities:pðΦÞ ¼

j pð/ Þ; pð/ Þ ¼ N ð/ j0; IÞ;

ðA1Þ

j pðu Þ; pðu Þ ¼ N ðu j0; IÞ;

ðA2Þ

n

n

n

n

pðUÞ ¼nnnnpðVÞ ¼j pðu Þ; pðu Þ ¼ N ðu j0; Λ Þ;
n

n

n

S

ðA3Þ

n

pðΛS Þ ¼ const:

and hyperparameters χ and β are used for the mixing matrices tohelp learn their dimension:pðCÞ ¼j pðc Þ; pðc Þ ¼ N c j0; ðk Þ v
kl

kl

kl



ðA4Þ

S k l

kl

pðBÞ ¼j pðb Þ; pðb Þ ¼ N b j0; ðk Þ b
km

km

km



S k m

ðA5Þ

km

Computation of the current model above is intractable due tothe joint probability of the parameters and factors. The variational
approximation is used, which restricts the joint posterior to a product of factor distributions, but allows the solution to be computed
analytically. The VB-EM algorithm iteratively maximizes the freeenergy F with respect to (w.r.t.) each factorized distribution to, atleast, a local maximum of F, alternating w.r.t. the posteriors q(U|Y)
and q(B|Y). Therefore, the following variational approximations
are made to make the model computationally tractable:
pðU; BjYÞcqðU; BjYÞ ¼ qðUjYÞqðBjYÞ
pðΦ; U; CjYÞcqðΦ; U; CjYÞ ¼ qðΦ; UjYÞqðCjYÞ

ðA6Þ

The update rules for the two-stage procedure are given. The
following posterior estimates are obtained for the factors in thefirst-stage VBE-step:qðUjYÞ ¼



n

n

n

P



j qðu jy Þ; qðu jy Þ ¼ N u j u ;Γ
1 P T

un ¼Γ

B ΛS yn ;

n

n

n

PT

P

P

n

Γ ¼ B ΛS BþKΨ1
B þI

ðA7Þ

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

In the first-stage VBM-step, the full posterior distribution of the
background mixing matrix B is computed, including its precision
matrix ΨB, and the MAP estimates of the noise precision Λ S andthe hyperparameter β.
qðBjYÞ ¼

j qðb jYÞ; qðb jYÞ ¼ N b j b ;ðk Þ Ψ
k

k



P

k

k

S k

B

k

P

B ¼RYU ΨB ; ΨB ¼ ðRUU þ bÞ1
1 PT Pb1 ¼ diag
B ΛS BþΨBK
1Pdiag RYY  BRTYU Þ
Λ1
S ¼
N

ð

Þ

ðA8Þ

Now that B and Λ S have been learned from the data, thestatistics of these noise sources are assumed not to change.
The second-stage VBE-step results as: P  P
PTqðkn jyn Þ ¼ N kn j kn ;Γ ; kn ¼Γ1Ω ΛS yn ;
P

PT

P

ðA10Þ

P



k

k

k

P

k

S k



C

P

1 ΨC ¼ ðRΦΦ þ cÞ
C ¼ðRYΦ  BRUΦ ÞΨ
C ;
1
P
P
T
c1 ¼ diag
C ΛS CþΨC
K

ð

þ

Þ

ðA11Þ

Solving for Λ̂−1 1, and further simplifying to:

RΦU ¼

P PT

fn un þN ΣΦU

RΦΦ ¼

n¼1

P PT

fn fn þN ΣΦΦ

ðA12Þ

n¼1

N
AL
A N
1X
¼
logjΛ1 j ðY  FGfn ÞT Λ1 ðY  FGfn Þ
AG AG 2
2 n¼1



ΣΦΦ
ΣUΦ

ΣΦU
ΣUU

!
¼0

Using Eq. (B4) and defining S as

1 
RYY  RYΦ R1ΦΦ R ΦY þ Σ0 ;
N þm

ðB6Þ

then Ĝ can be written as:1
Gr ¼ FT S1 F FT S1 RYΦ R1ΦΦ

where Σ = Γ− 1 is specified as:
Σ¼

ðB4Þ

Now Λ−1 1 is a function of G (since G was assumed known
when taking the derivative above); an expression for Λ−1 1 not dependent on G is needed. To solve for G (assuming Λ−1 1 is known),take the derivative of L :

S¼
N
X

1
ðRYY  FGRΦY þ Σ0 Þ
N þm

ðB5Þ

The matrices, such as RUΦ, represent the posterior covariance
between the two subscripts
N
X

Þ

N
m
1
logjΛ1 j þ logjΛ1 j  TrðΣ0 Λ1 Þ ¼ 0
2
2
2
ðB3Þ

1

j qðc jYÞ; qðc jYÞ ¼ N c j c ;ðk Þ Ψ
k

ð

N
AL
A
1X
¼

ðyn  FGfn ÞT Λ1 ðyn  FGfn ÞAΛ1 AΛ1
2 n¼1

Λ̂1 ¼

In the second-stage VBM-step, the posterior distribution of theinterference mixing matrix C is updated, including its precision
ΨC, as well as the MAP value of the hyperparameter χ. Thus, theposterior distribution for C is:qðCjYÞ ¼

where Σ0 and ν are hyperparameters. The graphical model in Fig. 1
indicates that G r and Λ1r are independent, and we give a flat prioron G r.
We choose ν = K + 2 for the distribution to be normalizable.
Whereas Σ0 could be inferred by directly measuring the sample
covariance, instead VBFA is used on the pre-stimulus data (like inthe first stage of SEFA, but applied to the post-stimulus data). From
VBFA on the post-stimulus data, Λ0 is the diagonal sensor precision
and B0 is the interference mixing matrix, so Σ0 = (B0B0T + Λ−0 1)− 1.
To solve for Λ1 (assuming G is known), take the derivative ofthe likelihood:ðA9Þ

Γ ¼ ΩT ΛS ΩþI ¼ Ω ΛSΩþKΨ1 þ I
 PT 
 1

PP
ΨC 0
¼ C
þI
PT ΛS ðC BÞþK
0
0
B

939

ðB7Þ


ðA13Þ

Since the expression for G is now known, this value can beplugged into Eq. (B4) to find Λ1 independent of G.

Appendix B. Derivation of NSEFALoc1 estimates

Appendix C. Derivation of NSEFALoc2 estimates

For each scanned voxel, we consider the likelihood function
over all the known data and hidden parameters:

In the VBE-step of NSEFALoc2, p(x n|yn) is found by finding
the q(x n|y n) that maximizes the free energy F and therefore best
approximates p(x n|y n), where Θ = {G, Λ2, α}. The variational
approximation, similar to the SEFA model in Eq. (A6), that theparameters and variables are conditionally independent given thedata, is used.





 
Lr ¼ log p Y; Gr ; Λr1 ~log p YjΛr1 ; Gr þ log p Λr1N X
log p YjΛr1 ; Gr ¼log p yn jΛr1 ; Gr

ðB1Þ

n¼1





p yn jΛr1 ; Gr ¼ N yn jFGr fn ; Λr11 r Tr Σ0 Λ1
 r
 r

r m=2 2
p Λ1 ¼ W Λ1 jm; Σ0 ~jΛ1 j e

ZF ðq; ΘÞ ¼

ðB2ÞdX dA qðXjY; ΘÞqðAjY; ΘÞ½log pðYjX; A; ΘÞþ log pðXÞ þ log pðAjΛ2 ; αÞ log qðXjY; ΘÞ  log qðAjY; ΘÞ
ðC1Þ

940

J.M. Zumer et al. / NeuroImage 41 (2008) 924–940

Maximizing for q(X|Y) yields
log qðXjYÞ ¼ EqðAjYÞ ð log pðY; X; AjΘÞÞðC2Þ

It can be shown that q(X|Y) is also a Gaussian distribution. The
mean of a Gaussian is the value that makes the derivative zero andthe variance of the Gaussian is the slope of gradient, yielding: P pðxn jyn Þ ¼ N xn j xn ;Γ
PTPP
xn ¼Γ1 A Λ2 ðyn  FGfn ÞðC3Þ

Γ ¼A Λ2 A þ KΨ þ I
TTo find the MAP estimate of G, the derivative of the free energy
is taken:AF
A¼
EqðXjYÞ EqðAjyÞ ½log pðYjX; A; GÞ ¼ 0
AG AG
"

#
N
A
1X
ðyn  FGfn  Axn ÞT Λ2 ðyn  FGfn  Axn Þ ¼ 0EqðXjYÞ EqðAjYÞ 
AG
2 n¼1

ðC4Þ
Plugging in the value for Ā, we obtain:


1
1
P 
G ¼ FT Λ2 F FT Λ2 RYΦ  RYX Ψ1 RXΦ RΦΦ  RΦX Ψ1 RXΦðC5Þ
The rest of the updates are given in NSEFALoc2.

References
Attias, H., 1999. Inferring parameters and structure of latent variable models
by variational bayes. Proc. 15th Conf. Uncert. Art. Intell., pp. 21–30.
Baryshnikov, B.V., Van Veen, B.D., Wakai, R.T., 2004. Maximum likelihood dipole fitting in spatially colored noise. Neurol. Clin. Neurophysiol. 2004, 53.
Bunch, P., Hamilton, J., Sanderson, G., Simmons, A., 1978. A free response
approach to the measurement and characterization of radiographicobserver performance. J. App. Photo. 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
for coherent source region suppression. IEEE Trans. Biomed. Eng. 53,
1357–1363.
Dalal, S.S., Guggisberg, A.G., Edwards, E., Sekihara, K., Findlay, A.M.,
Canolty, R.T., Knight, R.T., Barbaro, N.M., Kirsch, H.E., Nagarajan, S.S.,
2007. Spatial localization of cortical time-frequency dynamics. Conf. Proc.
IEEE Eng. Med Biol. Soc. 1, 4941–4944.
Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R.M., 2004. Mapping
human brain function with MEG and EEG: methods and validation.
NeuroImage 23 (Suppl 1), S289–S299.
Daunizeau, J., Friston, K.J., 2007. A mesostate-space model for EEG and
MEG. NeuroImage 38, 67–81.
Dogandzic, A., Nehorai, A., 2000. Estimating evoked dipole responses in
unknown spatially correlated noise with EEG/MEG arrays. IEEE Trans.
Sig. Proc., 13–25.
Friston, K., Harrison, L., Daunizeau, J., Kiebel, S., Phillips, C., TrujilloBarreto, N., Henson, R., Flandin, G., Mattout, J., 2008. Multiple sparsepriors for the M/EEG inverse problem. NeuroImage 39, 1104–1120.
Gelman, A., Rubin, D.B., 1996. Markov chain Monte Carlo methods inbiostatistics. Stat. Meth. 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.
Ikeda, S., Toyama, K., 2000. Independent component analysis for noisy data—
MEG data analysis. Neural Networks 13, 1063–1074.
Jun, S.C., George, J.S., Paré-Blagoev, J., Plis, S.M., Ranken, D.M., Schmidt,
D.M., Wood, C.C., 2005. Spatiotem-poral Bayesian inference dipole
analysis for MEG neuroimaging data. NeuroImage 28, 84–98.
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.
Kiebel, S.J., Daunizeau, J., Phillips, C., Friston, K.J., 2008. Variational
Bayesian inversion of the equivalent current dipole model in EEG/MEG.
NeuroImage 39, 728–741.
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.
Nagarajan, S.S., Attias, H.T., Hild, K.E., Sekihara, K., 2007a. A probabilistic algorithm for robust interference suppression in bioelectromagnetic sensor data. Stat. Med. 26, 3886–3910.
Nagarajan, S.S., Attias, H.T., Hild, K.E., Sekihara, K., 2007b. A
probabilistic algorithm for robust interference suppression in bioelectromagnetic sensor data. Stat. Med. 26, 3886–3910.
Nummenmaa, A., Auranen, T., Hämäläinen, M.S., Jääskeläinen, I.P., Lampinen,
J., Sams, M., Vehtari, A., 2007a. Hierarchical Bayesian estimates of
distributed MEG sources: theoretical aspects and comparison of variationaland MCMC methods. NeuroImage 35, 669–685.
Nummenmaa, A., Auranen, T., Hämäläinen, M.S., Jääskeläinen, I.P., Sams,
M., Vehtari, A., Lampinen, J., 2007b. Automatic relevance determination
based hierarchical Bayesian MEG inversion in practice. NeuroImage 37,
876–889.
Pascual-Marqui, R.D., 2002. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Meth. Find. Exp.
Clin. Pharmacol. 24, 5–12 Suppl D.
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.
Sekihara, K., Nagarajan, S., Poeppel, D., Marantz, A., 2002a. Performanceof an MEG adaptive-beamformer technique in the presence of correlated
neural activities: effects on signal intensity and time-course estimates.
IEEE Trans. Biomed. Eng. 49, 1534–1546.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., Miyashita, Y., 2002b.
Application of an MEG eigenspace beamformer to reconstructing spatiotemporal activities of neural sources. Hum. Brain Mapp. 15, 199–215.
Snodgrass, J.G., Corwin, J., 1988. Pragmatics of measuring recognition memory:
applications to dementia and amnesia. J. Exp. Psychol. Gen. 117, 34–50.
Trujillo-Barreto, N.J., Aubert-Vázquez, E., Penny, W.D., 2008. Bayesian M/EEG
source reconstruction with spatio-temporal priors. NeuroImage 39, 318–335.
Vrba, J., Robinson, S.E., McCubbin, J., 2004. How many channels areneeded for MEG? Neurol. Clin. Neurophysiol. 2004, 99.
Zumer, J., Attias, H., Sekihara, K., Nagarajan, S., 2006. April Two probabilistic algorithms for MEG/EEG source reconstruction, IEEE ISBI:
From Nano to Macro.
Zumer, J.M., Attias, H.T., Sekihara, K., Nagarajan, S.S., 2007. A probabilistic algorithm integrating source localization and noise suppression
for MEG and EEG data. NeuroImage 37, 102–115.