NeuroImage 183 (2018) 698–715

Contents lists available at ScienceDirect

NeuroImage
journal homepage: www.elsevier.com/locate/neuroimage

Hierarchical multiscale Bayesian algorithm for robust MEG/EEG
source reconstruction
Chang Cai a, Kensuke Sekihara b, c, Srikantan S. Nagarajan a, *abcDepartment of Radiology and Biomedical Imaging, University of California, San Francisco, CA, 94143-0628, United States
Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45 Yushima, Bunkyo-ku, Tokyo, 113-8519, Japan
Signal Analysis Inc., Hachioji, Tokyo, JapanA R T I C L E I N F O

A B S T R A C T

Keywords:
Brain mapping
Magnetoencephalography
Electroencephalography
Bayesian

In this paper, we present a novel hierarchical multiscale Bayesian algorithm for electromagnetic brain imaging
using magnetoencephalography (MEG) and electroencephalography (EEG). In particular, we present a solution to
the source reconstruction problem for sources that vary in spatial extent. We deﬁne sensor data measurements
using a generative probabilistic graphical model that is hierarchical across spatial scales of brain regions andvoxels. We then derive a novel Bayesian algorithm for probabilistic inference with this graphical model. Thisalgorithm enables robust reconstruction of sources that have different spatial extent, from spatially contiguous
clusters of dipoles to isolated dipolar sources. We compare the new algorithm with several representative
benchmarks on both simulated and real brain activities. The source locations and the correct estimation of source
time courses used for the simulated data are chosen to test the performance on challenging source conﬁgurations.
In simulations, performance of the novel algorithm shows superiority to several existing benchmark algorithms.
We also demonstrate that the new algorithm is more robust to correlated brain activity present in real MEG andEEG data and is able to resolve distinct and functionally relevant brain areas with real MEG and EEG datasets.

1. Introduction
Mapping of the entire brain's activity in humans is an important undertaking in cognitive neuroscience research that seeks to understandneural mechanisms of complex human behaviors. It also has clinicalapplications in patients with brain tumors and epilepsy, where functional
brain mapping is useful to guide neurosurgical planning, navigation, andresection.
Two techniques currently exist for non-invasive brain mapping ofelectrophysiological activity in humans: electroencephalography (EEG)and magnetoencephalography (MEG). MEG and EEG are complementary
techniques that measure, respectively, the magnetic ﬁeld outside thehead and the scalp electric potentials produced by electrical activity inneural cell assemblies. Since they directly measure electrical brain activity from neural ensembles, these methods have superior temporal
resolution compared to PET or fMRI, thereby enabling studies of the
dynamics of neural ensembles that occur at typical time scales on theorder of tens of milliseconds.
To estimate brain source activity from EEG or MEG data, sourcereconstruction algorithms are necessary, which consists of solving aforward problem and an inverse problem. The forward problem computes the scalp potentials and external magnetic ﬁelds for a speciﬁc set ofneural current sources for a given sensor conﬁguration, brain anatomy,
head geometries, and volume conductor properties. The inverse problem
estimates the parameters of neural sources from MEG and EEG sensor
data and makes use of the forward problem computations. The estimationof spatial locations and timing of brain sources is still a challenging
problem because it involves solving for unknown brain activity acrossthousands of voxels from the recordings of just a few hundred sensors. In
general, there are no unique solutions to the inverse problem because
there are many source conﬁgurations that could produce sensor data that
can account for the sensor observations. This nonuniqueness is referred
to as the ill-posed nature of the inverse problem. Besides handling the illposed nature of EEG or MEG imaging, the inverse algorithms have toaddress the challenge of searching for true source signals while minimizing the many sources of noise that interfere with the true signals.
Electrical, thermal and biological noise as well as background room
interference can be present.
To overcome these challenges, researchers have proposed many
efﬁcient inverse problem algorithms which can broadly be classiﬁed into

* Corresponding author.
E-mail address: sri@ucsf.edu (S.S. Nagarajan).
https://doi.org/10.1016/j.neuroimage.2018.07.056
Received 4 August 2017; Received in revised form 12 July 2018; Accepted 23 July 2018
Available online 27 July 20181053-8119/© 2018 Published by Elsevier Inc.

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 1. Graphical models for (a) BMN, (b) Champagne, (c) Tree_Champagne. Variables dependent on time are inside dotted box; Variables independent of time areoutside bottled box. Variables in circles are unknown and learned from the model, and Variables in squares are known. N is the number of voxels, si denotes the ith
voxel time course, si;j is the j-th voxel's time course in i-th region.

two categories: model-based parametric dipole ﬁtting and whole-brain
source imaging methods. Dipole ﬁtting methods assume that a small
set of current dipoles can adequately represent an unknown source distribution, which is a direct way to estimate source parameters and has
properties of high resolution but low accuracy. This is because solving for
dipole parameters requires nonlinear optimization over a highdimensional parameter space with solutions having great sensitivity toinitialization due to the high probability of being a local minima. This isespecially a signiﬁcant problem when multiple dipoles are considered.
Furthermore, estimating the number of dipoles remains an intractableproblem.
An alternative approach is whole-brain source imaging methods
which do not require prior knowledge of the number of sources and cangenerally avoid the non-linear search in the high dimensional parameter
space (Sekihara and Nagarajan, 2008, 2015; Wu et al., 2016; Wipf andNagarajan, 2009). These methods apply voxel discretization over a wholebrain volume, and assume a source at each voxel and estimate the amplitudes (and orientation) of the sources by minimizing a cost function.
Imaging methods can be further classiﬁed into two classes: tomographic
reconstruction and spatial scanning techniques. Tomographic techniquesmodel the activity at all candidate source locations simultaneously.
Tomographic techniques include minimum-norm estimation (MNE)(H€
am€al€ainen and Ilmoniemi, 1984, 1994), dynamic statistical parametricmapping (dSPM) (Dale et al., 2000), and standardized low -resolution
brain electromagnetic tomography (sLORETA) (Pascual-Marqui, 2002).
Some tomographic techniques promote sparseness in the solution (Matsuura and Okabe, 1995; Uutela et al., 1999), where the majority of thecandidate locations do not have signiﬁcant activity (Wipf and Nagarajan,
2009; Zumer et al., 2007, 2008; Wipf et al., 2010a). Empirical evidenceshows that a sparse source model can improve the accuracy of thelocalization in a noisy environment (Wipf et al., 2010a). In contrast,
spatial scanning techniques sequentially estimate the time course at
every candidate location while suppressing the interference from activity
at the other candidate source locations. Some examples of scanningtechniques are minimum-variance adaptive beamforming (Robinson
et al., 1992; Spencer et al., 1992; Van Veen et al., 1997; Sekihara et al.,1996) and other variants of beamformers (Sekihara and Nagarajan,
2008).
Most of the source reconstruction algorithms from the above classes
can be viewed in Bayesian framework (Sekihara and Nagarajan, 2015).
This perspective is useful because at a high level, the prior distribution,
implicitly or explicitly imposed, can be used to differentiate and compare
the various source localization methods. Recently, we have developedChampagne, a novel tomographic source reconstruction algorithm that isderived in an empirical Bayesian and incorporates deep theoretical ideas
about sparse-source recovery from noisy, constrained measurements.
Champagne improves upon existing methods of source reconstruction interms of reconstruction accuracy, robustness, and computational efﬁciency (Wipf et al., 2010a). Experiments with preliminary simulated andreal data, presented in (Owen et al., 2012), show that compared to othercommonly-used source localization algorithms, Champagne is morerobust to correlated sources and noisy data. However, when faced withmore complex brain activity patterns that span multiple spatial scales,such as clusters of dipolar sources or mixtures of clusters and isolated
dipolar sources, there are still no efﬁcient source reconstruction
algorithms.
Here, we present a novel hierarchical multiscale generative model forelectromagnetic measurements such as MEG and EEG. This algorithm canbe considered as a hierarchical multiscale extension of the Champagne
algorithm. We ﬁrst assume that brain voxels cluster into either anatomically or functionally deﬁned brain regions or parcels with region-levelspeciﬁc variances. Voxel activity is then assumed to have a component
arising from regions with additional voxel -speciﬁc variances to accountfor variations in voxel activity within a region. The voxel activity is thenassumed to be related to sensor data using standard lead-ﬁeld kernels
that are known given the geometry of the sensor measurements and thevolume conductor model. We then derive Bayesian algorithm for estimating voxel and region variances from sensor data. We present a novel
algorithm with both voxel and region variances, referred to as tree_Champagne. We evaluate its performance in simulations and real-data
and compare with existing benchmark algorithms.
2. Methods
This section describes the tree_Champagne algorithm including theprobabilistic generative model, estimation of the source and region activity, learning of hyperparameters, and its relation to other Bayesianinference algorithms.
2.1. The probabilistic generative model
We assume that MEG/EEG data have been collected for evoked or
induced source activity paradigms, with separate time-windows for
evoked or induced source activity and for background brain activity
including interference from biological, environmental sources and sensor
noise.
The generative model for the sensor data is:yðtÞ ¼

NX
li si ðtÞ þ ε

(1)

i¼1

where, yðtÞ 2 ℝdy 1 , is the output data of sensors at time t, dy is thenumber of channels measured, N is the number of voxels under consideration and li 2 ℝdy dc is the lead-ﬁeld matrix for the i-th voxel. The k-th
column of li represents the signal vector that would be observed at thescalp given a unit current source/dipole at the i-th voxel with a ﬁxed
orientation in the k-th direction. It is common to assume dc ¼ 2 (for MEG)
or dc ¼ 3 (for EEG), which allows ﬂexible source orientations to be
estimated in 2D or 3D space. Multiple methods based on the physical699

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 2. Example of the localization results for simulated MEG data with 3 clusters at SNIR ¼ 10. The activity power is normalized by the lead-ﬁeld value at each voxel.
The ground truth is shown for comparison.

700

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 3. Example of the localization results for simulated MEG data with 3 regions active at 10 dB. The activity power is normalized by the lead-ﬁeld value at eachvoxel. The ground truth is shown for comparison.

si ðtÞ ¼ vi ðtÞ þ gj zj ðtÞ

properties of the brain and Maxwell's equations are available for thecomputation of each li (Hallez et al., 2007). And si ðtÞ 2 ℝdc 1 is the i-th
voxel intensity at time t, which we assume it with dc orientations. Finally,
ε is a noise-plus-interference term where we assume, for simplicity, thatthe columns are drawn independently from Nð0; Σε Þ with known
covariance Σε . Temporal correlations can easily be incorporated if
desired using a simple transformation outlined in (Friston et al., 2008a)
or using the spatio-temporal framework introduced in (Bolstad et al.,2009). Here, we assume that the noise covariance can be estimated fromthe baseline and evoked data using a Stimulus-Evoked Factor Analysis,
SEFA (Nagarajan et al., 2007a) or variational Bayesian factor analysis(VBFA) model (Nagarajan et al., 2007b).
In our hierarchical framework, we divide the brain into R apriori
regions (or tiles) speciﬁed either anatomically or functionally (Tzourio-Mazoyer et al., 2002). The j-th region contains pj voxels. As a ﬁrst step,
we assume that the division of regions are assumed to be
non-overlapping, where each voxel belongs to exactly one region, but
this is not a necessity in the framework. Regional tiling may correspondto a map of anatomical or functional areas, or be constructed by, e.g.,
dividing the voxels into regions centered at equally-spaced locations
throughout the brain (Attias, 2013). We also assume that each unknown
region's activity zj ðtÞ 2 ℝdc 1 at time t is an equivalent dc -dimensional
neural current dipole, projecting from the j-th region. We then assume
that a given voxel's activity arises from the addition of the region's activity and voxel activity that is independent of the region's, as shown
below.

(2)

In the equation above, vi expresses the component that is intrinsic to
the i-th voxel and independent from activities of other voxels or the region a voxel belongs to. gj is the gain matrix between j-th region distribution and voxel si , here we assume it to be p1j , where pj is the number ofvoxels for j-th region. Then, the source data model in Eq. (1) is expressedsuch thatyðtÞ ¼
where

!
N
R
N
R
X
X
X
X
1X
li vi ðtÞ þli zj ðtÞ þ ε ¼
li vi ðtÞ þ
lj zj ðtÞ þ ε
pj i2ωj
i¼1
j¼1
i¼1
j¼1
P
i2ωj

(3)

indicates the summation regarding the voxels that belong to

the j-th region. We then denote the mean lead ﬁeld over the j-th region byP
lj : lj ¼ p1j
li and deﬁne an extended (voxel-augmented) lead ﬁeld matrixi2ωj

H such that
H ¼ l1 ; …; lN ; l1 ; …; lR ¼ ½h1 ; …; hNþR (4)

where hi ¼ li for i ¼ 1; …; N and hi ¼ liN for i ¼ N þ 1;…;N þ R. We alsodeﬁne an extended voxel vector, such that

T 
T
xðtÞ ¼ vT1 ðtÞ; …; vTN ðtÞ; zT1 ðtÞ; …; zTR ðtÞ ¼ xT1 ðtÞ; …; xTNþR ðtÞ(5)where xi ðtÞ ¼ vi ðtÞ for i ¼ 1; …; N and xi ðtÞ ¼ ziN ðtÞ for i ¼ N þ 1; …;
701

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 4. Example of the Aggregate Performance metric calculation with increasing correlation in clusters from 0.1 to 1 at 10 dB for 50 simulations: (A) Averaged Hit
Rate for all algorithms; (B) Averaged False Rate for all algorithms; (C) Averaged correlations for all hit sources; (D) Averaged Aggregate Performance scores for
all algorithms.

N þ R. Eq. (3) can then be rewritten asyðtÞ ¼ HxðtÞ þ εpðXjYÞ ¼

k¼1

(6)

ϒ1
6 0
ϒ ¼6
4 ⋮
0

0
ϒ2
⋮
0

3
⋯
0
⋯
0 7
7
⋱
⋮ 5
⋯ ϒ NþR

xk ¼ Γ 1 H T Σ1
ε yk

(12)

1

(13)

where
Σy ¼ Σε þ Hϒ H T(14)This Σy is called the model data covariance matrix. The solution in Eq.
(13) can be expressed as

(7)

3 2
ϒ1
x1 ðtk Þ
6 x2 ðtk Þ 7 6 0
7¼6
6
5 4 ⋮
4
⋮
0
xNþR ðtk Þ
2

(8)

0
ϒ2
⋮
0

32 T 3
⋯
0
h1
6 T 7
⋯
0 7
76 h2 7Σ1 y k⋱⋮ 54 ⋮ 5 y⋯ ϒ NþR
hTNþR

(15)

We can then express the source activity in terms of a spatial ﬁlter asshown below:

Using the noise assumption that ε  N ðεj0; Σε Þ, the conditional
probability pðYjXÞ is expressed as

k¼1

(11)

xk ¼ ϒH T ðΣε þ Hϒ H T Þ yk ¼ ϒH T Σ1
y yk

k¼1

K
K
Y
Y
pðYjXÞ ¼
pðyk jxk Þ ¼
N ðyk jHxk ; Σε ÞΓ 1 ¼ ϒ 1 þ H T Σ1
ε H

The posterior mean can be written in an alternative way, such that

The prior distribution is expressed asKYN ðxk j0; ϒÞpðXjϒÞ ¼

(10)

k¼1

where the variance and the mean are obtained as

The equation above is the data model used for the derivation of the
new algorithm.
The data vector yðtk Þ is denoted y k and the extended voxel vectorxðtk Þ is denoted xk for simplicity, tk is the time point at k. We formulate
the source reconstruction problem as the spatio-temporal reconstruction,i.e., the voxel time series x1 ; x2 ; …; xK is reconstructed using the sensor
time series y1 ; y2 ; …; y K . We express the whole time series x1 ; x2 ; …; xK
collectively as X 2 ℝðNþRÞdc K , and the whole time series y1 ; y 2 ; …; y K as
Y 2 ℝdy K .
We then deﬁne ϒ i as a prior variance dc  dc matrix of xi and deﬁne ϒ
as dc ðN þ RÞ  dc ðN þ RÞ block diagonal matrix expressed as
2

K
K
Y
Y
 

pðxk jyk Þ ¼
N xk xk ; Γ 1

xj ðtk Þ ¼ ϒ j hTj Σ1
y yk

(9)

(16)

2.3. Learning the hyperparameters ϒk¼1Here, the noise covariance Σε can be estimated using SEFA (Nagarajan
et al., 2007a) or VBFA (Nagarajan et al., 2007b) and is assumed to be
known for simplicity and subsequent considerations.

The Bayesian estimate of xk is given as the voxel posterior mean in Eq.
(12) or (16). In order to compute xk in Eq. (16), we need to know thehyperparameter ϒ. The hyperparameter ϒ is obtained by maximizingpðYjϒÞ which is called the marginal likelihood (Wipf et al., 2010b). The
marginal likelihood pðYjϒÞ is expressed as follows (details of the derivation of this function can be found in Appendix A).

2.2. Estimation of the source and region activity
To estimate the source distribution X, we ﬁrst derive the posterior
distribution pðXjYÞ, which is given by702

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 5. Simulation results of Aggregate Performance with four different conﬁgurations at 10 dB and 0 dB: (A) and (B) show results for increasing dipoles time courses
correlation from the same cluster; (C) and (D) show results for increasing correlation between clusters; (E) and (F) show results for increasing the number of clusters;(G) and (H) show results for variations in the sizes of the clusters.

"
K
1X
ðyk  Hxk ÞT Σ1log pðYjϒÞ ¼ 
ε ðy k  Hxk Þ
K k¼1
#
N þR
X
T
1
þ
xj ðtk Þϒj xj ðtk Þ  logjΣy jF ðϒ; ΛÞ ¼ 
(17)

j¼1

j¼1

Although the optimum value of the hyperparameter ϒ is obtained by
maximizing log pðYjϒÞ, maximizing the right-hand side of the equation
 
above is difﬁcult due to the inclusion of the last term logΣy .
 
Since logΣy  is a concave function with respect to ϒ (Wipf and

(19)

j¼1

where,
log pðYjϒÞ  F ðϒ; ΛÞ(20)always holds, and increasing F ðϒ; ΛÞ with respect to ϒ should result inincreasing the marginal likelihood log pðYjϒÞ. Therefore, the updatevalue of ϒ is derived as

Nagarajan, 2009), using dc  dc auxiliary-parameter matrices, Λj ðj ¼ 1;…; N þ RÞ, the relationship (Jordan et al., 1999) (Boyd and Vandenberghe, 2004),
 
 
NþR


X
 
T
tr Λj ϒ j  Λ0  logΣy 
 
j¼1

"
K
1X
ðyk  Hxk ÞT Σ1ε ðy k  Hxk Þ
K k¼1
#
NþR
NþR


X
X
xTj ðtk Þϒ 1xðtÞtr ΛTj ϒ j þ Λ0þjkjb ¼ argmaxF ðϒ; ΛÞ
ϒ
ϒ

(18)

Update rules can then be derived usingholds where Λ0 is a scalar term. Accordingly, we deﬁne an auxiliary costfunction F ðϒ; ΛÞ such that

703

(21)

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 6. Aggregate Performance with three different conﬁgurations: (A) and (B) show results for increasing the brain's regions at 10 dB and 0 dB; (C) and (D) show theperformance of all algorithms with ﬁxed 5 dipoles while increasing the number of clusters at 10 dB and 0 dB; (E) and (F) show results with ﬁxed 5 clusters but
increasing the number of dipoles at 10 dB and 0 dB.

∂

∂ϒ j



F ϒj; Λ

"
¼ ϒ 1j#K
1X
Txj ðtk Þxj ðtk Þ ϒ 1
j þ Λj ¼ 0
K k¼1

Maximum (EM) based algorithms such as those in (Friston et al., 2008a;
Sato et al., 2004).

(22)

Setting the right-hand side to zero, we get the equation,"ϒ j Λj ϒ j ¼

2.4. Algorithm summary#K
1X
xj ðtk ÞxTj ðtk Þ
K k¼1

Tree_Champagne is a source reconstruction algorithm based on
generative model Eq. (1) and is able to combine sparsity (from voxel level
inference) and smoothness (from regional-level inference) to produce
optimally smooth and sparse solutions.
Using the updating rules above, we can calculate the variance of both
voxels and regions. We denote variance of the voxel intrinsic component
as ϒ V and variance of the regions as ϒ R , the relationship between ϒ andϒ V , ϒ R is(23)

A positive semi-deﬁnite matrix that satisﬁes Eq. (23), can be derived
such that
#
"
"
#1=2
K
1X
1=2
1=2
T
b j ¼ Λj1=2 Λ1=2
x
ðt
Þx
ðt
Þ
Λ
Λj
ϒ
j k
j
j
j k
K k¼1

(24)

Eq. (24) is the update rule for ϒ j .
The update rule for Λj is derived using a fact that the hyperplanePNþRbT
j¼1 trð Λ j ϒ j Þ  Λ0 forms a tightest upper bound of the concave function logΣy . Such a hyperplane is found as the plane that is tangential to logΣy . Therefore, the update equation for Λj is derived as  b j ¼ ∂ logΣy  ¼ hT Σ1 hj
Λ
j y
 
∂ϒ

0
ϒR

(26)

3
ϒ1 ⋯ 04ϒ ¼ ⋮ ⋱ ⋮ 5
2 0 ⋯ ϒN3
ϒ Nþ1 ⋯
0R
ϒ ¼4 ⋮⋱⋮ 50⋯ ϒ NþR

(27)

ϒ¼

ϒV
0

where
2

V

(25)

j

In summary, the hyperparameter ϒ j are estimated by iterating Eq.
(16), Eq. (24) and Eq. (25). Each iteration is theoretically guaranteed toincrease (or leave unchanged) the cost function F ϒj ; Λ . The periteration cost is linear in the number of N þ R so the computational
cost is relatively modest (it is quadratic in dy , and cubic in dc , but thesequantities are relatively small). The convergence rate is orders ofmagnitude faster (Wipf and Nagarajan, 2009) than Expectation

The variance of i-th voxel is treated as the summation of the varianceof the i-th voxel's intrinsic component and variance of the region where
the i-th voxel belongs to. The time course of tree_Champagne can be
expressed as:704

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 7. Example of the localization results for 2 clusters and 2 dipoles at SNIR ¼ 10 dB. The activity power is normalized by the lead-ﬁeld value at each voxel. The
ground truth is shown for comparison.

705

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 8. EEG simulation results of the A Prime Metric (left column) and Aggregate Performance (right column) with three different conﬁgurations at 10 dB: (A) and (B)show results for increasing number of clusters; (C) and (D) show results with ﬁxed 5 clusters and increasing the number of dipoles; (E) and (F) show results with ﬁxed 5
dipoles while increasing the number of clusters.

2.5. Algorithm initialization
Initialization of the parameter updates are described here. First, Σε islearned from the pre-stimulus period using SEFA (Nagarajan et al.,2007a) or VBFA (Nagarajan et al., 2007b) and ﬁxed. Initialization for ϒ isset by running Bayesian Minimum-Norm (BMN) (Sekihara and Nagarajan, 2015), described below, to determine a whole-brain level varianceparameter and the variance of all voxels and regions are initialized to thisvalue. The precision and the mean of the posterior distribution pðxjyÞ arecomputed using Eq. (11) and Eq. (16). The hyperparameter ϒ is updated
using Eq. (24) and with the values of Λ updated using Eq. (25) and x
obtained earlier. Finally, we calculate the variance of both voxels andregions, the time course of each voxel using Eq. (28) and Eq. (29) withboth voxels and regions taken into consideration.
2.6. Relationship of tree_Champagne to BMN and to Champagne
BMN (Sekihara and Nagarajan, 2015) and Champagne are two other
Bayesian algorithms for source reconstruction which have close relation
to tree_Champagne. The difference among the three algorithms are in thegenerative model, as can be seen in Fig. 1. For BMN, voxels in source
space have a scalar variance ν. Bayesian estimation of this model yields
the BMN algorithm which results in smooth widespread activity
throughout the brain. In contrast to BMN, each voxel in Champagne has adifferent prior variance. Bayesian inference of the Champagne modelyields very sparse reconstructions (Wipf et al., 2010a). In contrast to bothof these algorithms, tree_Champagne uses a source space that issegmented into different regions according to prior anatomy and function. We then assume that each region has its own region's level variance.
Additionally, tree_Champagne also includes an intrinsic variance for eachFig. 9. Averaged Radius of Clusters with EEG simulations for Champagne,
tree_Champagne. The Ground Truth is shown for comparison.
tree T 1
bs treei ðtk Þ ¼ ϒ i hi Σy y k(28)

is expressed aswhere i-th voxel's variance ϒ treei
ϒ tree¼ ϒ Vi þ ϒ Rji(29)

where the i-th voxel belongs to j-th region.

706

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 10. Sensory Evoked Field localization results. The activity power is normalized by the lead-ﬁeld value at each voxel. All six algorithms localize to somatosensory
cortical areas, where Champagne and tree_Champagne are the most focal. BMN_sLORETA also performs well on the localization. Here we set the threshold fortree_Champagne and Champagne much lower than other benchmarks.

This section describes the performance evaluation of tree_Champagne
under different speciﬁc complex conﬁgurations compared with otherfour representative benchmark source reconstruction algorithms. Then
we evaluate the performance of tree_Champagne using real MEG and EEG
datasets with different tasks.

magnitude. For real MEG data tests, we use BMN (Sekihara and Nagarajan, 2015) for learning the regularization parameter which is then usedin conjunction with sLORETA - we call this algorithm of using sLORETA
for normalization after BMN as BMN_sLORETA. We found slight improvements in real data using BMN_sLORETA when compared to sLORETA with a ﬁxed regularization that we used in our simulations (seeFig. 10). For Champagne, Beamformer and sLORETA, we use Matlab ﬁles(nuts_Champagne.m, nuts_LCMV_Vector_Beamformer.m and nuts_sLORETA.m) from NUTMEG (Dalal et al., 2004). For MSP, we use the exact
implementation of MSP as included in standard settings in SPM12(spm_eeg_invert.m).

3.1. Benchmark source localization algorithms

3.2. Quantifying performance

Four representative source localization algorithms we chose tocompare with the performance of tree_Champagne are: (1) an adaptive
spatial ﬁltering method, linearly constrained minimum variance beamformer (referred to as Beamformer) (Robinson et al., 1992; Spencer et al.,
1992; Van Veen et al., 1997; Sekihara et al., 1996), (2) a non-adaptive
weighted minimum-norm method, standardized low-resolution tomographic analysis (referred to as sLORETA) (Dale et al., 2000; Pascual-Marqui, 2002), and two Bayesian based algorithms - (3) Champagne(Wipf et al., 2010a) and (4) MSP (Friston et al., 2008b). In simulations,for sLORETA we ﬁx the regularization to be 1e-6 times the maximum
eigenvalue of the composite lead-ﬁeld. We did not ﬁnd much variation in
performance when we changed the regularization by 1-2 orders of

In order to evaluate the performance on simulated results, two features are quantiﬁed: localization accuracy and time course estimation
accuracy. We ﬁrst examine whether sources are correctly localized, thenmeasure if the source time courses are accurately reconstructed for thosesource locations. The occurrence of both hits rate and false positives aretaken into account for the evaluation. The free-response ROC (FROC)
curve is used as it allows for multiple hits and false positives in a single
0
image (Darvas et al., 2004). The A metric (Snodgrass and Corwin, 1988)
estimates the area under the FROC curve for one hit rate (HR ) and false
positive rate (FR ) pair, or in our case, for each simulation. If the area
under the FROC curve is large, then the hit rate is higher compared to thefalse positive rate.

voxel independent of the regional variance. This variance partitioning
enables it to produce source reconstructions with varying spatial extents,
as we will show below.
3. Performance evaluation on simulation and real data707

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 11. Auditory Evoked Field results for three subjects. The activity power is normalized by the lead-ﬁeld value at each voxel. The results from both Champagne andtree_Champagne are shown in the last two columns, which outperform the other benchmark algorithms shown in the ﬁrst to three columns.

Fig. 12. Audio-Visual data localization results from tree_Champagne. The activity power is normalized by the lead-ﬁeld value at each voxel. Tree_Champagne is able
to localize a bilateral auditory response at 100 ms after the simultaneous presentation of tones and a visual stimulus. For bilateral auditory activity, the results oflocations and time courses are shown in (A), (B). Tree_Champagne can localize an early visual response at 150 ms after the simultaneous presentation of tones andvisual stimulus shown in (C) and (D).
0A ¼

HR  F R 1
þ2
2

two metrics that capture both the accuracy of the location and timecourses of the algorithms into a single metric called the Aggregate Per0
formance (AP) (Owen et al., 2012). It combines the A , R, and HR using
the following equation:(30)

where HR is calculated by dividing the number of hits by the true numberof seeded sources and FR is calculated by dividing the number of false
positive by the maximum number of false positives for each algorithm.
0
Eq. (30) is a simple way to compute A metric in our prior paper (Sekihara, 2016). The correlation coefﬁcient between the seed and estimatedsource time courses for each hit is used to determine the accuracy of thetime courses. We then average the correlation coefﬁcients for all the hits
for each simulation run, which is denoted as R. Finally, we combine these

AP ¼


1 0
A þ HR R
2(31)

The HR is used as a weight for R since we only compute the correlation
coefﬁcient for the sources that are correctly localized. AP ranges from0 to 1, with higher numbers reﬂecting better performance.

708

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 13. Face-processing task (MEG) localization results for tree_Champagne. The activity power is normalized by the lead-ﬁeld value at each voxel. Tree_Champagne
can localize an early visual response around 100 ms after the presentation of a face stimulus, results with time courses shown in subplot (A). A later visual response
around 200 ms after the presentation of a face stimulus are shown in subplot (B). The novel algorithm can localize the bilateral activation in fusiform gyrus that isthought to be in FFA, shown in (C) and (D). The peak for the brain activity is around 170 ms after the presentation of a face stimulus, and the time courses are shown
next to brain activity ﬁgures in subplots (C) and (D).

We expand conﬁgurations with correlations between clusters (inter-clusters, βinter ) which deﬁne the voxel time courses correlation fromdifferent cluster and correlations in clusters (intra-cluster, βintra ) which
deﬁne the voxel time courses correlation among the same cluster. We
also tested the effect of activity with both clusters and point sources. In
summary, the conﬁgurations we tested are as follows:

3.3. MEG simulations
In this paper, we generate data by simulating dipole sources withﬁxed orientation. Damped sinusoidal time courses are created as voxel
source time activity and we then project the voxel activity to the sensors
using the forward model generated lead ﬁeld matrix. The lead ﬁeld isconstructed within the brain volume assuming a single-shell sphericalmodel (Hallez et al., 2007) as implemented in SPM12 (http://www.ﬁl.
ion.ucl.ac.uk/spm) at the default spatial resolution of 8196 voxels at
approximately 5 mm spacing. The time course is then partitioned intopre- and post-stimulus periods. In the pre-stimulus period (180 samples) there is only noise plus interfering brain activity, while in the
post-stimulus period (300 samples) there is also source activities ofinterest on top of statistical similarity distributed noise plus interferingbrain activity. The noise plus interfering activity consists of actual
resting-state sensor recordings collected from a human subject presumed to have only spontaneous brain activity and sensor noise. The
voxel level activity is then projected to the sensors through the leadﬁeld and the noise/interference is then added to achieve a desired
signal to noise ratio. The simulated data has 271 sensor recordings. The
locations for the active sources are chosen so that there is some minimum distance between sources (at least 15 mm) and a minimum distance from the center of the head (at least 35 mm) (Julia P Owen et al.,2012).
We could adjust both the signal-to-noise-plus-interference ratio(SNIR) and the correlations between the different voxel time courses(inter-dipole αinter ) to examine the algorithm performance on unknown
correlated sources and ﬁxed orientation. In this paper, SNIR and correlation between sources are deﬁned in the same way as is shown in ourprior work (Julia P Owen et al., 2012).
Similar to our prior work, we picked difﬁcult conﬁgurations that wehave tested for Champagne (Julia P Owen et al., 2012). Additionally, inthis paper, we extend our tests to sources with extended spatial extent,
i.e. source clusters and regions with more complex conﬁgurations(Edwards et al., 2010). A voxel source is a point dipolar source and acluster source is deﬁned as sources with several closely located dipolar
sources. A region is set a priori using atlases by dividing the whole brain
into regions deﬁned either anatomically or functionally, such as theAutomated Anatomical Labeling (AAL) (Tzourio-Mazoyer et al., 2002).

1 Correlation within cluster - We examine the ability to reconstructclusters with increasing correlation of sources from the same cluster.
We seed 5 clusters with 20 sources for each cluster. The correlation ofsources from the same cluster is set as 0.1, 0.3, 0.5, 0.7 and 0.9 - in thissituation, we set the correlation between clusters as 0.25.
2 Correlation between clusters - We examine the inﬂuence of correlation between clusters for the novel algorithm. We seeded 5 clusters
with 20 sources per cluster. We set the correlation between clusters as
0.1, 0.3, 0.5, 0.7 and 0.9 - the source time courses within each cluster
is set to have an intra-cluster correlation coefﬁcient of 0.5.
3 Number of clusters - We test the ability to localize distributed clusters
by simulating different numbers of clusters. We seed 1, 4, 7, 10, 13, 16
clusters with 20 sources for each cluster. These clusters correspond to20, 80, 140, 200, 260 and 320 voxels having nonzero activity. The
placement of the cluster center is seeded randomly and cluster consists of sources seeded within the 19 nearest neighboring voxels.
4 Effect of clusters' size - We assess the robustness to localize distributed
sources with different cluster sizes. We seed 5 clusters with 10, 16, 22,
28, 34 and 40 active dipoles per cluster, which correspond to 50, 80,
110, 140, 170 and 200 active voxels.
5 Number of regions - Since our novel algorithm is based on the distribution of voxels into regions, we also test the inﬂuence of different
sizes of the region divisions. Here, we set the number of regions as 8,
9, 32, 95, 108, 116, 285 and 291 to evaluate performance of thetree_Champagne algorithm. For these simulations, we ﬁx the activity
as arising from 5 clusters with 20 sources for each cluster.
6 Mixed conditions (clusters and sources) - We extend the previous
cluster analysis experiments to investigate the effect of having bothcluster and dipole activity. We choose to set the number of clusters
from 1, 4, 7, 11, 14 to 17 with additional activity from 5 dipoles.
Subsequently, we set the number of clusters as 5 and vary the numberof dipoles from 1, 4, 7, 10, 13, 16 to 19.
709

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 14. Results for face processing (EEG) from novel algorithm and benchmarks. The ﬁrst row is the average power mapping from 0 ms to 400 ms, the second andthird rows are for peak power activity at 100 ms and 170 ms separately. Thresholds are 1% of the maximum activation of the image for tree_Chamapgne and 10% of
the maximum activation of the image for benchmarks. The activity power is normalized by the lead-ﬁeld value at each voxel.

plot these averaged results with standard error bars. We show the plots ofmean AP, our Aggregate Performance metric. We also show examples of
localization results from single simulations, which complement ouraggregate results.

If not indicated otherwise, each of the experiments is conductedwith the following settings: the source time courses within each cluster
have an intra-cluster correlation coefﬁcient of βintra ¼ 0:5 and an interdipole correlation coefﬁcient of βinter ¼ 0:25. We make the correlations
within the clusters higher than between clusters because nearby voxels
are more plausibly correlated than voxels at a distance. For clusters,
we are both interested in whether a cluster is localized and whether
the extent of a cluster is accurately reconstructed. To assess the00localization of clusters, we use the A metric. The A metric is calculated for clusters by testing if there is a local peak within the known
extent of the cluster. To assess the accuracy of the extent of clusters,
we calculate the fraction of seeded voxels with power in or above 10th
percentile of all voxels. At the same time, the power of localized peaks
should be at least 0.1 percentile of the maximum power.
The results obtained using simulated data are averaged over 50
simulations for each of six conﬁgurations with SNIR ¼ 0 or 10 dB and we

3.4. EEG simulations
We also test the novel algorithm on simulated EEG data using a scalar
lead-ﬁeld computed for a three-shell spherical model in SPM12 (http://
www.ﬁl.ion.ucl.ac.uk/spm) at the default resolution resulting in 8196
voxels at approximately 5 mm spacing. The simulated EEG data has 120
sensor recordings. With this lead-ﬁeld, EEG data is simulated in the same
way as the MEG data, as described above. We repeat the detection ofmultiple clusters and mixed conditions (clusters and sources) experiments for EEG simulations.

710

C. Cai et al.

NeuroImage 183 (2018) 698–715

Fig. 15. Epilepsy Spikes results for 7 subjects. The results of best time point dipole ﬁtting are shown in the left-most column, the results of benchmarks are shown fromsecond to forth columns, the novel algorithm's results are shown in the last column. The activity power is normalized by the lead-ﬁeld value at each voxel.

BC, Canada) with 1200 Hz sampling rate. The lead ﬁeld for each subject
is calculated in NUTMEG (Dalal et al., 2004) using a single-sphere headmodel (two spherical orientation lead ﬁelds) and an 8 mm voxel grid.
Each column is normalized to have a norm of unity. The data is digitallyﬁltered from 1 to 160 Hz to remove artifacts and the DC offset is3.5. Real datasets
All the MEG data here was acquired in the Biomagnetic Imaging
Laboratory at University of California, San Francisco (UCSF) with a CTF
Omega 2000 whole-head MEG system from VSM MedTech (Coquitlam,711

C. Cai et al.

NeuroImage 183 (2018) 698–715

benchmark algorithms perform somewhat similarly for these simulations, but their performance is not as good as tree_Chamagne andChampagne. Nevertheless, increasing the correlation in each cluster also
improves the performance of all algorithms.

removed.
We ran tree_Champagne and all of the benchmark algorithms on ﬁvereal MEG data sets: 1. Somatosensory Evoked Fields (SEF); 2. Auditory
Evoked Fields (AEF); 3. Audio-Visual Evoked ﬁelds; 4. Face-processing
task; 5. Interictal spike data from patients with epilepsy spikes. The
ﬁrst four datasets have been reported in our prior publications using theChampagne algorithm, and details about these datasets can be found in(Wipf et al., 2010a; Owen et al., 2012). Novel data included in this paper
are interictal spikes from seven patients with epilepsy. These spikes wereidentiﬁed by trained MEG technologists in the Biomagnetic Imaging
Laboratory, and the peak time-point was localized using dipole ﬁtting
method. For Champagne and tree_Champagne, we choose a pre-spike
window from 350 ms to 250 ms as a baseline control period and thepost-spikes window is from 50 ms to 50 ms where the spikes time is at
0 ms time point.
The EEG data (128-channel ActiveTwo system) was downloaded fromthe SPM website (http://www.ﬁl.ion.ucl.ac.uk/spm/data/mmfaces) andthe lead ﬁeld was calculated in SPM8 using a three-shell spherical modelat the coarse resolution. The EEG data paradigm involves randomized
presentation of at least 86 faces and 86 scrambled faces, here we subtract
the averaged scrambled-faces data to the averaged faces data to study thedifferential response to faces versus scrambled faces (Henson et al.,2010), and the power is plotted on a 3-D brain. The EEG data has beenreported in our prior publication using the Champagne algorithm, anddetails about our analyses of this dataset can be found in (Owen et al.,2012).

4.1.2. Inﬂuence of the correlation between clusters
The second row of Fig. 5 shows the inﬂuence of increasing the correlation between clusters on algorithm performances. Increasing thecorrelation between clusters has little inﬂuence on the performance of allalgorithms at both 10 dB and 0 dB. Based on the AP metric, it is clear thattree_Champagne outperforms all benchmarks. Although Champagne is
not as good as tree_Champagne, it is the best among all benchmarks whencompared to Beamformer, sLORETA and MSP.
4.1.3. Inﬂuence of the number of clusters
In the third row of Fig. 5, we plot the number of clusters versus AP
metric at SNIR levels of 10 dB and 0 dB. All algorithms have the same
trend at both 10 dB and 0 dB, with the increasing number of clusters, theAP score decreases. Again, tree_Champagne outperforms all benchmark
algorithms. Champagne is not as good as tree_Champagne but better thanothers. For benchmarks, at 10 dB, sLORETA shows higher AP score thanBeamformer and MSP. While at 0 dB, all benchmarks perform at a similarlevel when the number of clusters is more than 4.
4.1.4. Effect of clusters' size
The results of all methods at both 10 dB and 0 dB in response toincreasing clusters' size are presented in the last row of Fig. 5. Performances of all algorithms do not show much change when the clusters'
size increases. From the AP plot, tree_Champagne outperforms all
benchmarks. Again, Champagne is very close to tree_Champagne withsuperior performance when compared to Beamformer, sLORETA andMSP.

4. Results
4.1. MEG simulations
Fig. 2 shows a representative example of localization results for a
MEG simulation with 3 clusters at SNIR ¼ 10 dB, compared with theground truth. Champagne can ﬁnd all three clusters but it estimates activity that is more focal than the true spatial extent of the sources.
Tree_Champagne is also able to localize three clusters with estimates thatare more spatially distributed than Champagne. Beamformer is unable toﬁnd the three clusters correctly. sLORETA can ﬁnd all three clusters
correctly but produces blurred and diffuse solutions. In contrast, MSP canﬁnd all three clusters but reconstructions are smoother than ground truth
and also estimates additional sources that are not present in thesimulations.
The performance for a second special case where 3 regions of themodel are speciﬁed to be active is shown in Fig. 3. Only tree_Champagne
is able to reconstruct the correct active region, showing the extended
activity corresponding to each region. Champagne localizes the active
region but treats the regions' activity as if they are arising from several
point sources. In contrast, sLORETA, Beamformer and MSP do notaccurately estimate the regions' active and show blurred and inaccurate
reconstructions.
Fig. 4 shows an example of the steps that go into the aggregate performance metric calculation. With the increase of correlation in clusters,
we ﬁrst calculate Hit Rate (subplot A) and False Rate (subplot B) using
the method from our prior work (Owen et al., 2012). Then, the correlation between hit sources and seeded time series is obtained as shown in
0
subplot C. At last, we calculate the A matric and Aggregate performanceusing Eq. (30) and (31). Aggregate performance across 50 simulations foreach of 6 conﬁgurations is reported. For subsequent performance evaluation ﬁgures we only show the AP metric.

4.1.5. Effects of increasing the number of regions in the generative model
The ﬁrst row of Fig. 6 shows the inﬂuence on localization methods byincreasing the number of regions. The whole source space is segmented
into different size regions and tested at 10 dB and 0 dB. Although thisshould only inﬂuence the performance of tree_Champagne algorithm, we
also show performance for the benchmarks for these speciﬁc simulation
data instantiations using the same performance metrics. The intra-cluster
correlation is at 0.5 and the inter-clusters correlation is 0.25. The results
are averaged over 50 simulations each with 5 clusters seeded with 20
sources for each cluster, and the error bars show the standard error. As
we can see in the AP metric, when increasing the number of the regions,despite some changes in the performance of tree_Champagne, it is superior to the benchmark algorithms.
4.1.6. Performance for mixed source conﬁgurations
A single representative simulation experiment with 2 clusters and 2
dipoles at 10 dB is presented in Fig. 7, where the ground truth is shownon the ﬁrst row for comparison. Champagne and tree_Champagne canlocalize all clusters and the dipoles. Beamformer, sLORETA and MSP can
localize almost all activities but with very diffuse reconstructions andsome false positive activity estimates.
In order to evaluate the performance of source localization algorithms
for conﬁguration with both clusters and dipoles, we ﬁrst ﬁx the numberof dipoles as 5 and increase the number of clusters, then we ﬁx thenumber of clusters as 5 and increase the number of dipoles. The ﬁnal
results are plotted in the second and third rows of Fig. 6. As we can see in
the AP value with ﬁxed number of dipoles and increase number ofclusters, the performance of all algorithms decreases. Tree_Champagne
outperforms the benchmarks at both 10 dB and 0 dB. Champagne shows
better performance than other benchmarks at 10 dB but is close to others
at 0 dB. We then ﬁx the number of clusters and increase the number ofdipoles, and the performance of all algorithms decline as the number ofdipoles increases. Tree_Champagne still produces the highest scores

4.1.1. Inﬂuence of the correlation within each cluster
The sensitivity to performance as a result of increasing the correlation
within each cluster on both 10 dB and 0 dB is presented in the ﬁrst row ofFig. 5. From the AP plot, for both SNIR ¼ 10 dB and SNIR ¼ 0 dB, tree_Champagne outperforms all benchmarks. Champagne is not as good as
tree_Champagne but is much better than other benchmarks. The
712

C. Cai et al.

NeuroImage 183 (2018) 698–715

associated with the auditory stimulus. Tree_Champagne is able to localizebilateral auditory activity in Heschl's gyrus in the window around theM100 peak, shown in the ﬁrst row of Fig. 12. The two auditory sources
have the maximum power in the window around the M100 peak. We
show the early visual response in the second row of Fig. 12. Tree_Champagne is able to localize a source in the medial, occipital gyrus witha peak around 150 ms. We plot the power in the window around this peak
and the time course of the source marked with the cross hairs. Our novel
algorithm can localize a later visual response with a time course that has
power extending past 150 ms, which is similar to the results that we haveobtained with Champagne (Owen et al., 2012).

among all source localization algorithms both at 10 dB and 0 dB.
4.2. EEG simulations
In Fig. 8 we show EEG simulation results at 10 dB. According to ourtests, the performance of all algorithms have a similar trend with SNIR
equals to 10 dB or 0 dB. The left column shows the results of A Prime
Metric and the right column is the Aggregate Performance score. Acrossboth the A Prime and Aggregate Performance metrics, tree_Champagne
outperforms all benchmarks for all three different conﬁgurations.
In simulations, according to the evaluation function used in the paper,
the performance of tree_Champagne is much better than Champagne,
especially for clusters localization. Tree_Champagne is also more accurate than Champagne at estimating the spatial extent of cluster sources.
As is shown in Fig. 9, when we compare the radius of estimated size ofclusters for Champagne and tree_Champagne, the latter is better at estimating the spatial extent of the cluster.

4.4.4. Face-processing task: MEG
Localization of Face-processing task (MEG) in response to faces areshown in Fig. 13. We see an early visual cortical response to the presentation of the face visual stimulus in medial occipital cortex and later
visual cortical response more lateral to the early response shown in theﬁrst row of Fig. 13. Subsequently, tree_Champagne is able to localize thebilateral activation in the fusiform gyrus with peaks around 170 ms(Henson et al., 2010; Kanwisher et al., 1997). Performance of benchmarks algorithms on this dataset can be found in (Julia P Owen et al.,2012).

4.3. Summary for simulations
As we can see from the simulation results and analysis above, both at
10 dB or 0 dB, tree_Champagne outperforms all the benchmark sourcereconstruction algorithms. Next, we extend the evaluation of the performance using real MEG and EEG data.

4.4.5. Face-processing task: EEG
In Fig. 14, we present the results from using novel algorithm andbenchmarks on the face-processing task EEG data set. Fig. 14 shows theaverage power, M100 peak power and M170 peak power at different
rows separately. We see that tree_Champagne is able to localize the brain
activity with sparse peaks at visual areas and fusiform gyrus. However,the benchmarks produce the brain activity with either wrong location or
blurred solutions. Even though the threshold we use is 1% of the
maximum activation of the image for tree_Chamapgne and 10% of the
maximum activation of the image for benchmarks, our novel algorithm
gives us more sparse and accurate results.

4.4. Results of real data
This section shows the evaluation for our algorithms using real MEG
and EEG data, which contains ﬁve different MEG datasets and one EEG
dataset: Somatosensory Evoked Field Paradigm, Auditory Evoked Field,
Audio-Visual task, Face-processing task for MEG, Epileptic spikes data for
MEG and Face-Processing task for EEG.
4.4.1. Somatosensory Evoked Field Paradigm
Fig. 10 shows the results of the somatosensory evoked ﬁeld response
due to somatosensory stimuli presented to a subject's right index ﬁnger,
average derived from a total of 240 trials. A peak is typically seen 50 ms
after stimulation in the contralateral (in this case, the left) somatosensory
cortical area for the hand, i.e., dorsal region of the postcentral gyrus.
MSP, Champagne and tree_Champagne can localize this activation to thecorrect area of somatosensory cortex with focal reconstructions. Here, weshow performance in three benchmarks - Beamformer, sLORETA with aﬁxed regularization, and BMN_sLORETA. While benchmarks are also able
to localize somatosensory cortex, these reconstructions are more diffuse
especially for sLORETA with a ﬁxed regularization.

4.4.6. Epilepsy spikes
The localization results for epilepsy spikes data from seven patientsare shown in Fig. 15. The best time point dipole ﬁtting for each spike isshown in the left-most column for reference. As we can see, bothChampagne and tree_Champagne are able to localize almost all spikes forall subjects (shown in the forth to ﬁfth columns in Fig. 15). Champagne
and tree_Champagne perform similarly for all subjects. For otherbenchmark algorithms, Beamformer can localize the spike for each subject, but localization results are only reasonable for subject 5 since therest are either diffuse or have many spurious activations which arestronger than the true location of the spikes. BMN_sLORETA performs
better than Beamformer, but shows more diffuse results when compared
to Champagne and tree_Champagne. Since default MSP settings were
optimized for scalar lead-ﬁelds but these data included vector lead-ﬁelds,
we did not run MSP on these data.

4.4.2. Auditory Evoked Fields
The localization results for AEF data from three subjects are shown in
Fig. 11. The power of at each voxel in a 50–75 ms window around M100
peak is plotted for every algorithm. Both Champagne and tree_Champagne are able to consistently localize bilateral auditory activity for all
subjects (shown in the last two columns in Fig. 11). The activity is in
Heschl's gyrus, which is the location of primary auditory cortex. Champagne and tree_Champagne perform similarly for all subjects. Beamformer can ﬁnd the two auditory cortices only in one subject, whereas forthe rest of the subjects the activations are mostly biased towards thecentra of the head; This suggests that the correlation of bilateral auditory
cortical activity really impacts the performance of Beamformer.
BMN_sLORETA is able to ﬁnd the auditory activity for almost every
subject, but the results are diffuse and with additional spurious activities(not seen on the slices shown). MSP can localize bilateral auditory activity but with some location bias and more diffuse activation.

5. Discussion
This paper derives a novel hierarchical multiple spatial scale Bayesianalgorithm, tree_Champagne, for electromagnetic brain imaging using
magnetoencephalography (MEG) and electroencephalography (EEG)with comparisons to existing benchmark algorithms. The novel algorithm
is based on a principled cost function that maximizes the marginal likelihood of the data with fast, convergent update rules. The multiscale
formulation enables tree_Champagne to optimally combine smoothness(from regional-level inference) and sparsity (from voxel level inference).
Results show signiﬁcant theoretical and empirical advantages over many
existing methods. The algorithm readily handles multiple correlated
sources and is appropriate for sources that have variable spatial extent
ranging from isolated dipoles and extended clusters of dipoles, situations
that commonly arise even with simple cognitive neuroscience tasks.

4.4.3. Audio-visual evoked ﬁelds
Fig. 12 shows results of the audio-visual evoked ﬁelds for tree_Champagne. In subplot (A) and (B) we show the brain activations713

C. Cai et al.

NeuroImage 183 (2018) 698–715

using these update rules that are based on approximate likelihood
maximization using restricted maximum likelihood based, cost functions
that have slower convergence rates, tree_Champagne uses faster update
rules based on convex-bounds on true marginal likelihood of the data
(Wipf and Nagarajan, 2009). Given the similarities between the proposed
generative model, MSP and related hierarchical algorithms, inclusion of
proposed implementation ideas into these frameworks may minimize
observed differences in results.
Notably, however, the algorithms described in this paper do notincorporate temporal smoothness constraints and this represents the
future directions for our work. Various forms of temporal prior information or constraints can be uniﬁed within the framework of covariance component estimation. We are currently investigating the use of
temporal-smoothness priors in the form of basis functions (Zumer
et al., 2008) and in the form of autoregressive smoothness priors, which
also model spatiotemporal correlations in the background noise and canpotentially improve performance. The best example of such an effort isthe Bayesian Electromagnetic Spatio-Temporal Imaging of Extended
Sources (BESTIES) (Liu et al., 2016) algorithm, which is built upon a
Bayesian framework that determines the spatio-temporal smoothness ofsource activities in a fully data-driven fashion is based on a Markov
Random Field (MRF), which can precisely capture local cortical interactions, employed to characterize the spatial smoothness of source
activities, and importantly the temporal dynamics of which aremodeled by a set of temporal basis functions (TBFs). Jean Daunizeau
et al. (Jean et al., 2006) also introduced a Bayesian framework toincorporate distinct temporal and spatial constraints on the solution
and to estimate both parameters and hyperparameters of the model. A
full multivariate autoregressive (MAR) model formulates directed interactions (i.e., effective connectivity) between sources. The observation process of MEG data, the source dynamics, and a series of the priors
are combined into a Bayesian framework using a state-space representation. By formulating the source dynamics in the context of MEG
source reconstruction, and unifying the estimations of source amplitudes and interactions, the effective connectivity without requiring theselection of regions of interest can be identiﬁed (Fukushima et al.,2015). We derive inspiration for our future work from these approaches. Our next steps will focus on incorporating
temporal-smoothness and effective connectivity prior on our novel algorithms, which hold promise for improving upon an already robust
source localization algorithm.

The experiments with simulated data exemplify that tree_Champagne
provides robust localization and time course estimation with complex
source conﬁgurations and noisy data for both MEG and EEG simulations
with correlated sensor data. Tree_Champagne outperforms existing
benchmarks with highly correlated sources even at high levels of interference at 0 dB. We also found that with increasing of the number ofclusters and the size of clusters, tree_Champagne performs much better
than the benchmark algorithms. Notably, tree_Champagne performance
better than Champagne since tree_Champagne shows more extended
activity for clusters. For more complex conﬁgurations with simultaneousclusters and dipoles activity for both MEG and EEG, tree_Champagne isalso able to accurately localize the simulated activity and signiﬁcantly
outperforms benchmark algorithms.
Experiments with real data highlight the source localization abilitiesof the novel algorithm. It is difﬁcult to evaluate localization accuracy
with real data since the ground truth is not known. For this reason, we
have chosen real data sets that have well-established patterns of brain
activity; AEF, audio-visual, and face-processing data. For all these realdata, the tree_Champagne algorithm performs superiorly compared tobenchmarks and improves upon our prior work on Champagne. Additionally, here we examine a novel dataset of interictal spikes from patients with intractable epilepsy. For these data, tree_Champagne is able tosuccessfully localize all spikes for all subjects.
In this paper, the novel algorithm mainly models and addresses issues
related to incorporating priors for spatial-smoothness of sources activity.
We extend our prior framework of Champagne to include this spatialsmoothness using regional variances. Other researchers have taken
different approaches for incorporating priors on spatial-smoothness forsources reconstruction. Kn€
osche (Kn€osche et al., 2013) has proposed afunctional similarity as priors for the reconstruction of distributed source
current densities from EEG: patchLORETA1, which uses both topologicalneighborhood and prior information to deﬁne smoothness and patchLORETA2, which neglects topological neighborhood (Fukushima et al.,2015). Alternatively, fMRI-Informed Regional Estimation (FIRE) (Ouet al., 2010) utilizes information from fMRI in EEG/MEG source reconstruction which takes advantage of the spatial alignment between theneural and vascular activities, while allowing for substantial differencesin their dynamics.
The region-based variance model in tree_Champagne is different fromother multiscale or hierarchical approaches in several ways (Ou et al.,2010; Friston et al., 2008c; Jean et al., 2006; Babadi et al., 2014). First,
algorithms like the Multiple Sparse Priors (MSP) also evaluated here
(Friston et al., 2008a) impose spatial kernel smoothness across voxels
based on the adjacency matrix and only include regional level varianceswith no voxel-level variances. Second, in contrast to these algorithms, we
do not use variational approximations to factorize the posterior variancesat the region-level and voxel-levels, which allows for the posterior voxel
and regional variances to be correlated. Finally, we do not use greedy
algorithms like those proposed in Babadi et al. and Friston et al. (Friston
et al., 2008c; Babadi et al., 2014), which are highly sensitive to initialization and have the possibility of sub-optimal solutions. In contrast to

Acknowledgment
The authors would like to thank Susanne Honma, Danielle Mizuiri
and Anne Findlay for collecting much of the MEG data in the Biomagnetic
Imaging Laboratory. We would also like to thank Julia Owen for the
support of MEG data simulation and evaluation, Hagai Attias for inspiration and early discussions, and Inez Raharjo for editing. This work was
supported by NIH grants R01EB022717, R01DC013979, UCOP MRP-17454755, and a gift from Ricoh Company Ltd.

Appendix A. Derivation of the marginal likelihood function
Here, we derive the expression for the marginal likelihood function shown in Eq. (17). We make the use of the form (Sekihara and Nagarajan, 2015)(pp.244)log pðYjϒÞ ¼ EpðXjYÞ

log pðY; XjϒÞ¼
pðXjYÞ

Z
dXpðXjYÞlog

pðY; XjϒÞpðXjYÞ

(A.1)¼ EpðXjYÞ ½log pðYjXÞ þ EpðXjYÞ ½log pðXjϒÞ þ H ðpðXjYÞÞ
Substitution of equations Eq. (8), Eq. (9) and Eq. (10) into Eq. (A.1) results in the relationship714

C. Cai et al.

NeuroImage 183 (2018) 698–715

hXK

i
T 1
log pðYjϒÞ ¼ logjΣε j  EpðXjYÞðyHxk Þ Σε ðy k  Hx k Þ
k
k¼1
hXK
i
xT ϒxk  logjΓjþlogjϒj  EpðXjYÞ
k¼1 k
#
"
K
NþR
X
1X
T 1
T
1
¼
xj ðtk Þϒ j xj ðtk Þðyk  Hxk Þ Σε ðyk  Hxk Þ þ
K k¼1
j¼1
þlog

(A.2)

jΣε jjϒjjΓjusing equation (4.28) from book (Sekihara and Nagarajan, 2015) (pp.55), we get equation Eq. (17).

Liu, Ke, Yu, Zhu Liang, Wu, Wei, Gu, Zhenghui, Li, Yuanqing, Nagarajan, Srikantan, 2016.
Bayesian electromagnetic spatio-temporal imaging of extended sources with Markov
Random Field and temporal basis expansion. Neuroimage 139, 385–404.
Matsuura, K., Okabe, Y., June 1995. Selective minimum-norm solution of the biomagnetic
inverse problem. IEEE (Inst. Electr. Electron. Eng.) Trans. Biomed. Eng. 42 (6),608–615.
Nagarajan, Srikantan S., Attias, Hagai T., Hild, Kenneth E., Sekihara, Kensuke, 2007.
A probabilistic algorithm for robust interference suppression in bioelectromagnetic
sensor data. Stat. Med. 26 (21), 3886–3910.
Nagarajan, S.S., Attias, H., Hild 2nd, K.E., Sekihara, K., 2007. A probabilistic algorithm
for robust interference suppression in bioelectromagnetic sensor data. Stat. Med. 20,
3886–3910.
Ou, Wanmei, Nummenmaa, Aapo, Ahveninen, Jyrki, Belliveau, John W.,
H€am€al€ainen, Matti S., Golland, Polina, 2010. Multimodal functional imaging using
fMRI-informed regional EEG/MEG source estimation. Neuroimage 52 (1), 97–108.
Pascual-Marqui, R.D., 2002. Standardized low resolution brain electromagnetic
tomography (sLORETA): technical details. Methods Find. Exp. Clin. Pharmacol. 24,
5–12.
Robinson, S.E., Rose, D.F., 1992. Current source image estimation by spatially ﬁltered
MEG. In: Hoke, M., et al. (Eds.), Biomagnetism Clinical Aspects. Elsevier Science
Publishers, pp. 761–765.
Sato, Masa-aki, Yoshioka, Taku, Kajihara, Shigeki, Toyama, Keisuke, Goda, Naokazu,
Doya, Kenji, Kawato, Mitsuo, 2004. Hierarchical bayesian estimation for MEG inverse
problem. Neuroimage 23 (3), 806–826.
Sekihara, K., 2016. Computing resolution for neuromagnetic imaging systems. J Comput
Eng Inf Technol 5, 3. https://doi.org/10.4172/2324 9307:2.
Sekihara, K., Nagarajan, S.S., 2008. Adaptive Spatial Filters for Electromagnetic Brain
Imaging. Springer-Verlag, Berlin, Heidelber.
Sekihara, Kensuke, Nagarajan, Srikatan S., 2015. Electromagnetic Brain Imaging: aBayesian Perspective. Springer-Verlag, Berlin, Heidelber.
Sekihara, K., Scholz, B., 1996. Generalized Wiener estimation of three-dimensional
current distribution from biomagnetic measurements. In: Aine, C.J., et al. (Eds.),
Biomag 96: Proceedings of the Tenth International Conference on Biomagnetism.
Springer-Verlag, New York, pp. 338–341.
Snodgrass, Joan G., Corwin, June, 1988. Pragmatics of measuring recognition memory:
applications to dementia and amnesia. J. Exp. Psychol. Gen. 117 (1), 34.
Spencer, M.E., Leahy, R.M., Mosher, J.C., Lewis, P.S., November 1992. Adaptive ﬁlters formonitoring localized brain activity from surface potential time series. In: Conference
Record for 26th Annual Asilomer Conference on Signals, Systems, and Computers,
pp. 156–161.
Tzourio-Mazoyer, Nathalie, Landeau, Brigitte, Papathanassiou, Dimitri, Crivello, Fabrice,
Etard, Olivier, Delcroix, Nicolas, Mazoyer, Bernard, Joliot, Marc, 2002. Automatedanatomical labeling of activations in SPM using a macroscopic anatomicalparcellation of the MNI MRI single-subject brain. Neuroimage 15 (1), 273–289.
Uutela, Kimmo, H€am€al€ainen, M., Somersalo, Erkki, 1999. Visualization ofmagnetoencephalographic data using minimum current estimates. Neuroimage 10(2), 173–180.
Van Veen, B.D., Van Drongelen, W., Yuchtman, M., Suzuki, A., 1997. Localization of brain
electrical activity via linearly constrained minimum variance spatial ﬁltering. IEEE
Trans. Biomed. Eng. 44, 867–880.
Wipf, David, Nagarajan, Srikantan, 2009. A uniﬁed bayesian framework for MEG/EEG
source imaging. Neuroimage 44 (3), 947–966.
Wipf, David P., Owen, Julia P., Attias, Hagai T., Sekihara, Kensuke, Nagarajan, Srikantan
S., 2010. Robust bayesian estimation of the location, orientation, and time course ofmultiple correlated neural sources using MEG. Neuroimage 49 (1), 641–655.
Wipf, D.P., Owen, J.P., Attias, H.T., Sekihara, K., Nagarajan, S.S., 2010. Robust Bayesianestimation of the location, orientation, and time course of multiple correlated neuralsources using MEG. Neuroimage 49, 641–655.
Wu, Wei, Nagarajan, Srikantan, Chen, Zhe, 2016. Bayesian machine learning: EEG/MEG
signal processing measurements. IEEE Signal Process. Mag. 33 (1), 14–36.
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.
Zumer, Johanna M., Attias, Hagai T., Sekihara, Kensuke, Nagarajan, Srikantan S., 2008.
Probabilistic algorithms for MEG/EEG source reconstruction using temporal basis
functions learned from data. Neuroimage 41 (3), 924–940.

References
Attias, Hagai T., 2013. M/EEG imaging by learning mean norms in brain tiles. In:
Biomedical Imaging (ISBI), 2013 IEEE 10th International Symposium on. IEEE,
pp. 548–551.
Babadi, Behtash, Obregon-Henao, Gabriel, Lamus, Camilo, H€am€al€ainen, Matti S.,
Brown, Emery N., Purdon, Patrick L., 2014. A subspace pursuit-based iterative greedy
hierarchical solution to the neuromagnetic inverse problem. Neuroimage 87,427–443.
Bolstad, Andrew, Van Veen, Barry, Nowak, Robert, 2009. Space time event sparse
penalization for magneto-/electroencephalography. Neuroimage 46 (4), 1066–1081.
Boyd, Stephen, Vandenberghe, Lieven, 2004. Convex Optimization. Cambridge university
press.
Dalal, Sarang S., Zumer, J.M., Agrawal, V., Hild, K.E., Sekihara, Kensuke, Nagarajan, S.S.,
2004. NUTMEG: a neuromagnetic source reconstruction toolbox. Neurol. Clin.
Neurophysiol. Annu. Vol.: NCN 52, 2004.
Dale, A.M., Liu, A.K., Fischl, B.R., Buckner, R.L., Belliveau, J.W., Lewine, J.D., Halgren, E.,
2000. Dynamic statistical parametric mapping: combining fMRI and MEG for highresolution imaging of cortical activity. Neuron 26, 55–67.
Darvas, Felix, Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R.M., 2004. Mapping humanbrain function with MEG and EEG: methods and validation. Neuroimage 23,
S289–S299.
Edwards, Erik, Nagarajan, Srikantan S., Dalal, Sarang S., Canolty, Ryan T., E Kirsch, Heidi,
Barbaro, Nicholas M., Knight, Robert T., 2010. Spatiotemporal imaging of corticalactivation during verb generation and picture naming. Neuroimage 50 (1), 291–301.
Friston, Karl, Harrison, Lee, Jean, Daunizeau, Kiebel, Stefan, Phillips, Christophe, TrujilloBarreto, Nelson, Henson, Richard, Flandin, Guillaume, Mattout, Jrmie, 2008.
Multiple sparse priors for the M/EEG inverse problem. Neuroimage 39 (3),1104–1120.
Friston, Karl, Harrison, Lee, Jean, Daunizeau, Kiebel, Stefan, Phillips, Christophe, TrujilloBarreto, Nelson, Henson, Richard, Flandin, Guillaume, Mattout, Jeremie, 2008.
Multiple sparse priors for the M/EEG inverse problem. Neuroimage 39 (3),1104–1120.
Friston, Karl, Chu, Carlton, Mour~ao-Miranda, Janaina, Hulme, Oliver, Rees, Geraint,
Penny, Will, Ashburner, John, 2008. Bayesian decoding of brain images. Neuroimage
39 (1), 181–205.
Fukushima, Makoto, Yamashita, Okito, Kn€
osche, Thomas R., Sato, Masa-aki, 2015. MEG
source reconstruction based on identiﬁcation of directed source interactions onwhole-brain anatomical networks. Neuroimage 105, 408–427.
Hallez, Hans, Vanrumste, Bart, Grech, Roberta, Joseph, Muscat, De Clercq, Wim,
Vergult, Anneleen, D'Asseler, Yves, Camilleri, Kenneth P., Fabri, Simon G., Van
Huffel, Sabine, Lemahieu, Ignace, Nov 2007. Review on solving the forward problem
in EEG source analysis. J. NeuroEng. Rehabil. 4 (1), 46.
H€
am€
al€
ainen, M.S., Ilmoniemi, R.J., 1984. Interpreting Measured Magnetic ﬁelds of theBrain: Estimates of Current Distributions. Technical Report TKK-F-A559. Helsinki
University of Technology.
H€
am€
al€
ainen, M.S., Ilmoniemi, R.J., 1994. Interpreting magnetic ﬁelds of the brain:minimum norm estimates. Med. Biol. Eng. Comput. 32 (35–42).
Henson, Richard N., Flandin, Guillaume, Friston, Karl J., Mattout, Jeremie, 2010.
A Parametric Empirical Bayesian framework for fMRI-constrained MEG/EEG sourcereconstruction. Hum. Brain Mapp. 31 (10), 1512–1531.
Jean, Daunizeau, Mattout, Jeremie, Clonda, Diego, Goulard, Bernard, Benali, Habib,
Lina, J.-M., 2006. Bayesian spatio-temporal approach for EEG source reconstruction:
conciliating ECD and distributed models. IEEE (Inst. Electr. Electron. Eng.) Trans.
Biomed. Eng. 53 (3), 503–516.
Jordan, Michael I., Ghahramani, Zoubin, Jaakkola, Tommi S., Saul, Lawrence K., 1999.
An introduction to variational methods for graphical models. Mach. Learn. 37 (2),
183–233.
Julia P Owen, David P Wipf, Attias, Hagai T., Sekihara, Kensuke, Nagarajan, Srikantan S.,
2012. Performance evaluation of the champagne source reconstruction algorithm on
simulated and real M/EEG data. Neuroimage 60 (1), 305–323.
Kanwisher, Nancy, McDermott, Josh, Chun, Marvin M., 1997. The fusiform face area: amodule in human extrastriate cortex specialized for face perception. J. Neurosci. 17(11), 4302–4311.
Kn€osche, Thomas R., Gr€aser, Markus, Anwander, Alfred, 2013. Prior knowledge on cortex
organization in the reconstruction of source current densities from EEG. Neuroimage
67, 7–24.

715