1

Robust Empirical Bayesian Reconstruction ofDistributed Sources for Electromagnetic Brain
ImagingChang Cai, Mithun Diwakar, Dan Chen, Member, IEEE Kensuke Sekihara, Fellow, IEEEand Srikantan S. Nagarajan, Fellow, IEEE

Abstract—Electromagnetic brain imaging is the reconstructionof brain activity from non-invasive recordings of magnetic fields
and electric potentials. An enduring challenge in this imaging
modality is estimating the number, location, and time courseof sources, especially for the reconstruction of distributed brain
sources with complex spatial extent. Here, we introduce a novelrobust empirical Bayesian algorithm that enables better reconstruction of distributed brain source activity with two key ideas:
kernel smoothing and hyperparameter tiling. Since the proposed
algorithm builds upon many of the performance features of thesparse source reconstruction algorithm - Champagne and we refer to this algorithm as Smooth Champagne. Smooth Champagne
is robust to the effects of high levels of noise, interference andhighly correlated brain source activity. Simulations demonstrate
excellent performance of Smooth Champagne when compared
to benchmark algorithms in accurately determining the spatial
extent of distributed source activity. Smooth Champagne also
accurately reconstructs real MEG and EEG data.
Index Terms—Electromagnetic Brain Mapping, Bayesian Inference, Distributed brain activity, Inverse problem, Magnetoencephalography, Electroencephalography

I. I NTRODUCTIONONINVASIVE functional brain imaging has made atremendous impact in improving our understanding of
the human brain. Functional magnetic resonance imaging(fMRI) has been the predominant modality for imaging thefunctioning brain during the past 2 decades. However, fMRI
lacks the temporal resolution required to image the dynamic
and oscillatory spatiotemporal patterns associated with activities in the brain. Direct non-invasive measurements ofthese neuronal oscillatory activity in the millisecond timescale can be achieved with magnetoencephalography (MEG)and electroencephalography (EEG), this temporal resolution
is essential for epileptic form activity localization and for
imaging dynamic of brain networks subserving perception,action and cognition.
MEG/EEG sensor data only provides qualitative information about underlying brain activities. Analysis is typically
performed based on the qualitative analyses of experiencedN

Chang Cai, Mithun Diwakar and Srikantan S. Nagarajan arewith Department of Radiology and Biomedical Imaging, University ofCalifornia, San Francisco, CA 94143-0628 e-mail: sri@ucsf.edu (seehttps://radiology.ucsf.edu/research/labs/biomagnetic-imaging).
Dan Chen is with School of Computer, Wuhan University, China.
Kensuke Sekihara is with Department of Advanced Technology in
Medicine, Tokyo Medical and Dental University, 1-5-45 Yushima, Bunkyo-ku,
Tokyo 113-8519, Japan, and Signal Analysis Inc., Hachioji, Tokyo.

users regarding the sensitivity profile of the sensors. To extract
more precise information from the sensor data, it is essential toreconstruct actual brain activity (source space reconstruction)from the recorded sensor data. This process involves two steps:
the forward and inverse models. The forward model uses
source, volume conductor, and magnetic field measurement
models to calculate the lead field or gain matrix that describes
a linear relationship between sources and the measurements.
Inverse algorithms are then employed to estimate the parameters of neural sources from MEG and EEG sensor data. The illposed inverse problem for source localization with MEG/EEG
data involves estimating brain activity from noise sensor data
where the number of brain voxels (typically 3000 to 10000)
is much larger than the number of sensors (typically ∼ 300
sensors for MEG and ∼ 128 electrodes for EEG).
A wide variety of source localization algorithms exist for
estimating source activity to overcome the difficult inverse
problem. These algorithms may be roughly organized into
three classes: dipole fitting, spatial scanning, and tomography.
Dipole fitting can produce very sparse results but with two
major caveats. First, the process of nonlinear optimization
can result in poor performance when multiple sources arepresent due to local minima in the solution space. Second, foradequate performance, the number of active neuronal sourcesmust be known a priori. Scanning algorithms, also known
as spatial filters, are based on reconstructing activity at a
specific voxel of a discretized cortex. Spatial filter methods
do not need a priori knowledge of the number of the brain
sources. However, traditional spatial filter methods such as
beamformers perform poorly when different regional brain
activities are correlated in time. In contrast, tomographic
methods are estimating activity at all voxels simultaneously.
Some tomographic methods are able to promote sparseness inthe solution [1], [2], [3], where the majority of the candidate
locations do not have significant activity. Empirical evidenceshows that a sparse source model can improve the accuracyof localization in a noisy environment.
Most source reconstruction algorithms from the three categories can be framed in a Bayesian schema [4]. This perspective is useful because it demonstrates that various source localization methods are interrelated and, in fact, can be thought ofas implicity or explicity imposed prior distributions within a
Bayesian framework. Algorithms such as MNE [5], [6], dSPM
[7], sLORETA [8], MVAB (and other beamformers) [9], [10],[11], [12], [13] assume a known, fixed prior. Alternatively,

2

the parameters of the prior distribution (hyperparameters) canbe learned from the data for robust performance, commonlyreferred to as empirical Bayesian algorithms [4]. One suchrecent empirical Bayesian source reconstruction algorithm that
we have developed is referred to as Champagne [1]. Champagne outperforms several benchmark algorithms in a variety
of simulated source configurations, especially for sparse brain
source activity [14]. Since the development of the originalChampagne algorithm, we have proposed several extensionsof this framework. We have proposed hierarchical spatiotemporal extensions of Champagne to deal with sparse spatialand smooth temporal sources [15]. Recently, we proposeda hierarchical version of the Champagne algorithm, called
tree Champagne, for reconstruction of sources with mixed
spatial extents at a voxel and at a regional-level as specified
by anatomical or functional atlases [16].
Here, we develop a novel algorithm designed for accurate
reconstruction of source activity with distributed spatial extents. Since the algorithm builds upon the Champagne framework, we refer to it as Smooth Champagne. We propose two
key ideas to better handle distributed source activity: kernelsmoothing and hyperparameter tiling. First, we introduce akernel filter that converts a spatially smooth generative model
for brain sources into a sparse subspace. Second, we assume
that this sparse subspace is clustered into regions or tiles,
and we estimated variance hyperparameters for this model
at a resolution of regions [1]. Model learning of SmoothChampagne deploys robust empirical Bayesian Inference anduses a principled cost function which maximizes a convex
lower bound on the marginal likelihood of the data. Resulting
update rules are fast and convergent. In Section II, we derive
this new algorithm and in Sections III and IV we demonstrate
its performance in simulated and real MEG and EEG data incomparison to benchmarks, followed by a brief discussion in
Section V.
II. T HEORY
A. The probabilistic generative model
The generative model for the sensor data with distributed
brain activity is defined asy(t) =NP

Li si (t) + ε = LS(t) + ε,

(1)

i=1

where, y(t) = [y1 (t), . . . , yM (t)]T is the output of the sensors
at time t, M is the number of channels, N is the numberof voxels under consideration and Li ∈ RM ×q is the leadfield matrix for i-th voxel, q is the number of directions
for each voxel. The k-th column of Li represents the signal
vector that would be observed at the scalp given a unit current
source/dipole at the i-th vertex with a fixed orientation in thek-th direction. It is common to assume q = 2 (for MEG) or
q = 3 (for EEG), which allows flexible source orientations tobe estimated in 2D or 3D space. Multiple methods based on thephysical properties of the brain and Maxwell’s equations areavailable for the computation of each Li [17]. si (t) ∈ Rq×1 isthe brain activity for i-th voxel at time t with q directions.
The term ε is the noise-plus-interference whose statistics

are estimated from the baseline period. We assume that thecolumns are drawn independently from N (0, Σ ε ), where Σ εis the noise covariance which could be estimated using avariational Bayesian factor analysis (VBFA) model [18]. Temporal correlations can easily be incorporated if desired using asimple transformation outlined in [19]. And L = [L1 , . . . , LN ]and S(t) = [s1 (t)T , . . . , sN (t)T ]T . For simplicity, we definematrix Y = [y(1), . . . , y(T )] and S = [S(1), . . . , S(T )] as theentire sensor and source time series. T is the number of time
points.
In order to reconstruct the distributed brain activity robustly,
we assume that distributed brain source activity arises froma sparse set of clustered sources that (1) are spatially contiguous and (2) share a common time courses. We define
the sparse brain activity with clustered sources at time t as
Z(t) = [z1 (t)T , z2 (t)T , . . . , zN (t)T ]T , where zi (t) ∈ Rq×1 .
We also define a local spatial smoothing kernel matrix as
H = [h1 , h2 , . . . , hN ]T and hi = [hi1 , hi2 , . . . , hiN ]T is the
local spatial smoothing kernel for i-th voxel. S(t) = HZ(t)ensures spatial contiguity. The model of Eq. (1) can then be
expressed asy(t) = LHH−1 S(t) + ε = GZ(t) + ε,(2)where G = LH is the new modified leadfield matrix, G =[G1 , . . . , GN ] and Gi = Li hi .
There are many ways to define the local spatial smoothingkernel matrix hi . For mathematical convenience here, we
define hi , a smoothing kernel that includes flexible parameterization for different kernel properties such as width, pass-band
and transition-band characteristic, as follows:
j = i,

p 
 1
dij
(3)
hij =
√
j 6= i, exp −
m D
where, we denote the minimum distance between i-th and j-th
(j 6= i) voxel as D and the distance between i-th voxel andj-th voxel is denoted as dij . In the above formulation, theparameter m determines the width of the smoothing kernel,
uniformity in the pass-band, and the parameter p determines
the slope of the transition region of the kernel. For the rest of
this manuscript, unless otherwise specified, the default values
for p is 4 and for m is 2.
Additionally, we divide the whole brain voxels into R
regions (or tiles) specified either anatomically or functionallyand the r-th region contains Nr voxels. By default the numberof regions is assumed to be around one tenth of the total
number of voxels, regions are assumed to be non-overlappingin this paper. We assume that voxels within the same brain
region share a common variance hyperparameter. The priordistribution for sources can be written asp(Z|Ω) =T Y
R YY
t=1 r=1 n∈ζr

N (zn (t)|0, Ωr I),(4)

where Z = [Z(1), . . . , Z(T )] is the entire sparse time series,ζr is the set of voxels in r-th region and n ∈ ζr indicatesthat the n-th voxel is in r-th region. Ωr is the prior variancehyperparameter of r-th region, I is a q×q identity matrix. If we3

use the distribution of noise and interference as ε ∼ N (0, Σ ε ),the conditional probability p(Y|Z) can be expressed asp(Y|Z) =TQ

p(y(t)|Z(t)) =

t=1

T
Q
t=1

N (y(t)|GZ(t), Σ ε ).

(5)

Eqs. (2), (4) and (5) form the sparse model appropriate for
distributed sources, next, we will use Bayesian inference toderive the update rules of model parameters and hyperparameters. In the limit of when the kernel width or the tile/region
size is equal to the number of voxels in the lead-field, thegenerative model will result in the algorithm commonly known
as Bayesian minimum-norm algorithm. In such a limit, thealgorithm will fail to capture distributed source clusters thatare on the spatial scale of 5-35 contiguous voxels. Conversely,at the other limit of when the kernel width and the tile/region
size is equal to that of a single voxel, the generative model
will become similar to that for the Champagne algorithm, i.e.
the algorithm will become more suited to reconstruct sparsesource activity.
B. Bayesian inference
To estimate the sparse source activity Z, we use Bayesianinference and derive the posterior distribution p(Z|Y), which
is given byp(Z|Y) =T
Q

p(Z(t)|y(t)) =

t=1

T
Q
t=1

N (Z(t)|Z̄(t), Γ −1 ),(6)where the mean and the variance [1], [16], [20] are(7)

Γ −1 = Υ + GT Σ ε G.

(8)

G

T

Υ is a qN × qN diagonal matrix expressed asΥ1 0 ···0 0 Υ2 ···0 Υ = .
.
..  ,
.
..
..
 ..
. 
0
0 ··· ΥN(9)

Υ i = Ωr I if i-th voxel belongs to r-th region.
In order to compute Z̄(t) in Eq. (7), we need to know
the hyperparameter Ωr . The hyperparameter Ωr is obtainedby maximizing the cost function F(Ω, Φ), with an auxiliary
variable Φ, and is expressed as
T
T

1 P
y(t) − GZ(t) Σ −1y(t) − GZ(t)
ε
T t=1
R
T
P
P
1
−
zn (t)T (Ωr I)−1 zn (t)T r=1,n∈ζr t=1
RP−
tr(Φr Ωr I) + Φ0 .

The optimization problem of Φr is equivalent to finding
the hyperplane that forms a closest upper bound of log |Σ y |.
Such a hyperplane is found as the plane that is tangential tolog |Σ y |. Therefore, we set the derivative of log |Σ y | withrespect to Υ to zero, the update rule for Φr is obtained as
Φr =


1 X
tr GTn Σ−1
y Gn .
q

(12)

n∈ζr

The update rule for zn (t) can also be derived by setting thederivative of the new cost function F(Ω, Φ) with respect tozn (t) to zero. The source time course for the n-th voxel atr-th region is expressed aszn (t) = Ωr IGTn Σ−1y y(t).

(13)

The proposed algorithm then repeats the update steps
Eqs. (11), (12) and (13) until the cost function Eq. (10)converges.
Finally, the estimated sparse brain activity for n-th voxel is
Eq. (13) and the distributed brain activity issn (t) = hn zn (t) = hn Ωr IGTn Σ−1y y(t).

Σ −1
ε y(t),

Z̄(t) = Γ−1

We then set the derivative of the cost function F(Ω, Φ)with respect to Ω to zero, the regional resolution level hyperparameter Ωr for r-th region is given by [20]:v
u 1 PTutr zn (t)zn (t)T
1 t T Nr t=1,n∈ζr
.
(11)
Ωr =
q
Φr

(14)

A summary of the proposed algorithm is shown in Table 1.
Table 1 The Smooth Champagne Algorithm
1: Input sensor data y(t), leadfield matrix L.
2: Create the local smoothing kernel matrix H using Eq (3),calculate the new modified leadfield matrix using G = LH.
3: Set appropriate initial values to Ωr (r = 1, . . . , R), andΣε.
4: Compute Φr , and zn (t), (n = 1, . . . , T ) using Eq. (12)and Eq. (13).
repeat5: Update Ωr using Eq. (11).
6: Update Φr , and zn (t) using Eq. (12) and Eq. (13).
until The cost function Eq. (10) converges.
7: Derive the time course using Eq. (14), output s(t).

F =−

(10)

r=1

Eq. (10) is a convex lower bound of the logarithm of marginal
likelihood of the data p(Y|Ω) (details of the derivation of thisfunction can be found in Appendix A).

III. P ERFORMANCE EVALUATION ON SIMULATION AND
REAL DATAThis section describes algorithm performance evaluations
for different simulated complex source configurations andreal datasets. We also compared performance with benchmark
source localization algorithms.

4

A. Benchmark algorithms for comparison
We chose the following representative source localization
algorithms to compare with the performance of SmoothChampagne: (1) an adaptive spatial filtering method, linearly
constrained minimum variance beamformer (referred to as
Beamformer) [9], [10], [11], [12], two non-adaptive weighted
minimum-norm method, (2) standardized low-resolution tomographic analysis (referred to as sLORETA) [7], [8], (3)a variant of mixed-norm algorithm (MxNE) specially tailoredto handle multiple time-points and unconstrained source orientations [2], [3], and two Bayesian based algorithms: (4)Champagne [1] and (5) MSP [19]. Champagne, Beamformer,
and sLORETA, were implemented using NUTMEG libraries(nuts Champagne.m, nuts LCMV Vector Beamformer.m andnuts sLORETA.m) [21]. Standard implementation of MSP
from SPM12 was used (spm eeg invert.m).
B. Choice of tuning parameters
In Smooth Champagne there are two free parameters –namely the width of the smoothing kernel and the tile/region
size. All other parameters are estimated from data. In principle,
these two parameters can be estimated from data by computing
the model evidence as a function of kernel width and thetile/region size. For other algorithms like MxNE, sLORETA
and Beamformer the main free parameter is the regularization
constant. In our implementation of MxNE, we set the regularization parameter based on the sensor noise levels frombaseline data. For sLORETA, we use the default setting in
NUTMEG [21]. For Beamformer, we set the regularization
parameter to be 1e-3 times the maximum eigenvalue of thesensor data covariance.
C. Initialization
The initialization of the parameters for Smooth Champagne
are as follows. We learn the Σ ε from the baseline period using
Variational Bayesian Factor Analysis [18] and then assume it
is fixed when estimating other model parameters from actual
data. The local spatial smoothing kernel is set by Eq. (3) andafter that, the value of hij less than 0.1 is set as 0 to eliminate
the influence of voxels that are too far away. Initialization forΩ is set by first running Bayesian Minimum-Norm (BMN)[20] to determine a whole-brain level variance parameter, thenthe variance of region level is set by the average variance ofvoxels from the same region.
D. Quantifying performance
To evaluate the performance of Smooth Champagne in simulations, we use free-response receiver operator characteristics(FROC) which shows the probability for detection of a true
source in an image vs. the expected value of the number offalse positive detections per image [22], [16], [14], [23]. Based0on the FROC, we compute an A metric [24] which is an
estimate of the area under the FROC curve for each simulation.
If the area under the FROC curve is large, then the hit rate is0higher compared to the false positive rate. The A metric is
computed as follows:
0

A =

1
HR − FR
+
2
2

(15)

We make modifications to this performance metric to be applicable for reconstruction of distributed and clustered sourcesthat are balanced metrics across a variety of sparse andnon-sparse reconstruction algorithms. We use the following
definitions of hit rate and false rate for the performance
evaluation. First, the voxels localized by each algorithm thatare included in the calculation of hit-rates are defined as voxelsthat are i) at least 0.1% of the maximum activation of thelocalization result; and ii) within the largest 10% of all of thevoxels in the brain.
Within these subset of voxels, we test whether each voxel
is within the three nearest voxels to a true source cluster. If aparticular voxel’s estimated activity lies within a true cluster,that voxel gets labeled as a ’sub-hit’. If the number of ’sub-hits’
within a particular true cluster is more than 50% of the numberof voxels in that cluster, that cluster gets labeled as a ’hit’.
HR is then calculated by dividing the number of hit clusters
by the true number of clusters. For FR , the false positive
voxels localized by each algorithm are defined as voxels thatlie within the largest 10% of the whole voxels. FR is definedby dividing the number of potential false positive voxels by
the total number of inactive voxels for each simulation. Lastly,0
we use Eq (15) to calculate the extent scores A metric whichranges from 0 to 1, with higher numbers reflecting better
performance.
E. MEG simulations
In this paper, we generate data by simulating dipole sources
with fixed orientation. Damped sinusoidal time courses are created as voxel source time activity and then projected to the sensors using the leadfield matrix generated by the forward model,
assuming 271 sensors and a single-shell spherical model [17]as implemented in SPM12 (http://www.fil.ion.ucl.ac.uk/spm)
at the default spatial resolution of 8196 voxels at approximately 5 mm spacing and unless otherwise indicated, we set
the number of voxels in a region to be around eight resulting in
about one thousand regions, and the average kernel width is set
to around ten voxels. The time course is then partitioned intopre- and post-stimulus periods. In the pre-stimulus period (480samples) there is only noise and baseline brain activity, whilein the post-stimulus period (480 samples) there is also source
activities of interest on top of statistically similarity distributed
noise plus interfering brain activity. The noise activity, whichis added to achieve a desired signal to noise ratio, consists ofactual resting-state sensor recordings collected from a humansubject presumed to have only spontaneous brain activity andsensor noise.
Signal-to-noise ratio (SNR), correlations between voxel
time courses in the same cluster (intra-cluster αintra ), andcorrelations between voxel time courses from different clusters (inter-clusters αinter ) were varied to examine algorithm
performance. SNR and time course correlation are defined in a5standard fashion [16], [14]. The following configurations were
tested:1 Cluster number: 3, 5, 7, 11, 13, and 15 clusters wereseeded with 20 sources each, corresponding to 60, 100,
140, 220, 260, and 300 active voxels. Each cluster
consists of the center, chosen randomly, and 19 nearest
neighbors.
2 Cluster size: 5 clusters were seeded with 5, 10, 15, 20,
25, 30, 35 and 40 active dipoles each.
3 SNR: Simulations were performed with 5 randomly
seeded clusters at SNRs of 0 dB, 3 dB, 6 dB, 9 dB,
12 dB, and 15 dB.
4 Inter-cluster correlation: 5 clusters were seeded with 20
sources each. Correlation between clusters was varied
from 0.1, 0.3, 0.5, 0.7, and 0.9. Intra-cluster correlation
was set to 0.8.
5 Influence of width of local smoothing kernel and sizeof tiles: We evaluate the algorithm performance as afunction of local smoothing kernel width. We examine
performance for reconstruction of 5 clusters of different
sizes (from 5-35 voxels per cluster). We increase thewidth of the local smoothing kernel from 3 voxels to
121 voxels while the tile size is maintained at 6 voxels.
We also examine performance as a function of average
tile size ranging from 4-132 voxels per tile and fix theaverage kernel width (15 voxels).
If not indicated above, each of the simulations was conducted with the following parameters: SNR of 10 dB, intracluster correlation coefficient of αintra = 0.8, and inter-cluster
correlation coefficient of αinter = 0.25. Correlations withinclusters were modeled higher than between clusters to simulate0
real cortical activity. Results for A and localization wereobtained by averaging over 50 simulations for each of theconfigurations above.
F. EEG simulations
The algorithm was also tested on simulated EEG data witha scalar lead field computed for a three-shell spherical modelin SPM12 (http://www.fil.ion.ucl.ac.uk/spm) with more than
5000 voxels at approximately 8 mm spacing and 120 sensors.
EEG sensor recordings were simulated through the forward
model in the same way as the MEG data. Since the algorithm
has similar performance with EEG data as with MEG, only
one specific example is shown in the paper. Three clusters areseeded with SNR of 20 dB, intra-cluster correlation of 0.9,
and inter-cluster correlation of 0.9.
G. Real datasets
Real MEG data 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, BC, Canada) with 1200 Hz samplingrate. The lead field for each subject was calculated in NUTMEG [21] using a single-sphere head model (two sphericalorientation lead fields) and an 8 mm voxel grid. Each column
was normalized to have a norm of unity. The data weredigitally filtered from 1 to 160 Hz to remove artifacts andDC offset.
Six algorithms were run on three real MEG data sets: 1. Somatosensory Evoked Fields (SEF); 2. Auditory Evoked Fields(AEF); 3. Audio-Visual Evoked fields. The data sets havebeen reported in our prior publications using the Champagne
algorithm, and details about these datasets can be found in [1],
[14].
The
EEG
data
(128-channel
ActiveTwo
system) was downloaded from the SPM website(http://www.fil.ion.ucl.ac.uk/spm/data/mmfaces)andthelead field was calculated in SPM8 using a three-shell
spherical model at the coarse resolution of 5124 voxels
at approximately 8 mm spacing. The EEG data paradigm
involves randomized presentation of at least 86 faces and
86 scrambled faces. The averaged scrambled-faces data wassubtracted from the averaged faces data to study differential
response [25]. The EEG data has been reported in our priorpublication using the Champagne algorithm, and details about
our analyses of this dataset can be found in [16], [14].
IV. R ESULTS
A. Simulation Results
1) MEG simulations: Figure 1 shows a representative example of localization results for an MEG simulation with 3
clusters. In this configuration, the SNR, correlation of dipole
activities within the cluster and between clusters are 10 dB,
0.9 and 0.9, respectively. Compared with the ground truth,
Beamformer and sLORETA can localize all three clusters
but produce blurred and spurious activity, while MxNE canlocalize two of the three clusters but produces focal sources
and misses one cluster. In contrast, MSP can localize all
three clusters, but reconstructions are smoother than groundtruth and produces additional sources that are not present
in the simulations. Champagne can accurately localize all
three clusters but estimates activity that is more focal than
the true spatial extent of the sources. Smooth Champagne isable to localize three clusters which are the closest spatiallydistributed to the ground truth.
1 Influence of the number of clusters - In Figure 2 (A), we0plot the number of clusters versus A metrics at SNRlevels of 10 dB and correlations within and between
clusters as 0.8 and 0.25. All algorithms have the same
trend, with decreasing performance as number of clusters
increases. However, Smooth Champagne outperforms all0benchmark algorithms. MxNE shows higher A scoresthan MSP, Beamformer and sLORETA.
2 Effect of increasing cluster size - Results of all algorithms
in response to increasing cluster size are presented in
Figure 2 (B). From the performance scores, SmoothChampagne outperforms all benchmarks except for when
the size of clusters is less than 10 voxels. In this sparsesource activity scenario, as to be expected, Champagne
performs better than Smooth Champagne. Champagne
performs better than other benchmark algorithms whencluster size is less than 35, when cluster size increases
to 40, MSP outperforms Champagne. In contrast, MSP,Cluster 1

6

Ground Truth

Beamformer
Beamformer

sLORETA

Beamformer

sLORETA

MxNE

MSP

MxNE

MSP

MxNE

Champagne
Champagne

MSP

Smooth Champagne
6

Smooth Champagne

Champagne

Smooth Champagne

Cluster 1

Ground Truth

Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Sub-lobar (L)Insula (L)White Matter (L)Left Cerebrum
Parietal Lobe (L)Precuneus (L)White Matter (L)Left Cerebrum
Parietal Lobe (L)Precuneus (L)White Matter (L)Left Cerebrum
Parietal Lobe (L)Precuneus (L)White Matter (L)Left Cerebrum
Parietal Lobe (L)Posterior Cingulate(L)
White Matter (L)Left Cerebrum
Parietal Lobe (L)Precuneus (L)White Matter (L)Left Cerebrum
Parietal Lobe (L)Precuneus (L)White Matter (L)Left Cerebrum
Parietal Lobe (L)Precuneus (L)White Matter (L)Cluster 2
Cluster 3

Cluster Cluster
3
3

Cluster 3

Cluster 2

Cluster 2Cluster 2 Cluster 1 Cluster 1

Ground TruthsLORETA

1: A single MEG simulation example of the localization results with three clusters for six algorithms: Beamformer,
MCE,
MSP,
Champagne
SmoothThe groundtruthof clusters
shown
for algorithms:
comparison. Beamfo
Fig.
1: A
single MEG andsimulation
examplethe localization
results
withlocation
three clusters
for
six arealgorithms:
Beamformer,
AETA,
singleMEG
simulation
exampleofChampagne.
theof localization
resultswiththreeclusters
for six
SmoothFig.
Champagne
capturesthetrue spatial
extentof alltheThesources.
In thiscorrelationBeamformer,
ofsLORETA,MCE,
MSP, simulation
Champagne
andSmoothgroundtruth
locationof clusters
are thefor algorithms:
comparison.
1: A single
MEG
exampleofChampagne.
thelocalization
results
withconfiguration,
threeclusters
for SNR,
six
A,
MCE, MSP,
Champagne
andSmoothChampagne.
Thegroundtruth
locationof shownclusters
are shownfor compa
Only SmoothChampagne
captures thetrue spatial
extentof 0.9all thesources.
In this configuration, the SNR, correlation ofe activities
withintheclusterandbetweenclusters
are
10
dB,
and
0.9,
respectively.
sLORETA, MCE, MSP, Champagne and Smooth Champagne. The ground truth location of clusters are shown for comparison.
Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Right Cerebrum
Posterior Lobe (R)Declive (R)Right Cerebrum
Posterior Lobe (R)Declive (R)ooth Champagne
capturesthetrueandspatial
extentsources.
In this forconfiguration,
the SNR, correlati
dipole activities
within thecluster
between clusters
areof10 all
dB, the
0.9 and0.9, respectively.
Fig.
1: A single
MEG simulation
the localization
three clusters
six algorithms:
Only
SmoothChampagne
captures example
the trueofspatial
extent ofresults
all thewith
sources.
In this configuration,
theBeamformer,
SNR, correlation oftivitiesdipole
withinthe MxNE,
clusterandcluster
betweenclusters
are 10aredB,
0.9
andand0.9,
respectively.
sLORETA,MSP,
Smooth clusters
Champagne.
Thegroundlocationof clusters are shown for comparison.
activities
withintheChampagne
andandbetween
10 dB,
0.9truth0.9, respectively.
Only Smooth Champagne captures the true spatial extent of all the sources. In this configuration, the SNR, correlation of dipole
activities within the cluster and between clusters are 10 dB, 0.9 and 0.9, respectively.
0.9

0.9

0.85

0.85

0.8

0.8
0.75

0.9

0.9
0.7

A PrimeA Prime

0.75

0.85

0.9

0.7

0.9

0.85

demonstrate improved performance with increasing SNR.
Again, Smooth Champagne outperforms all benchmark
0.8
0
0.8
algorithms. MxNE shows higher A scores than MSP,
0.8
0.7
0.6
sLORETAand Beamformer.
cluster size
cluster size
0.6
0.7
5
5
4
Influenceofthecorrelation
betweenclusters
The
subplot
0.5
15
15
Fig.Simulation
3: Simulation
results ofof thethe performanceof the
novel
0.6
0.4
Fig.
3:
results
performance
of
the
novel
25
25
0.4
(D)
of
Figure2shows
theinfluenceofincreasing
thealgorithm
asafunctionoflocalsmoothingkernel
width
35
35and
0.5
algorithm ascorrelationa functionof localsmoothingkernelperformance.
width and
0.3
betweenclusters
onforalgorithm
0.2
the
tile
size.
We
examineperformancereconstruction
of
0.4
0.2
Fig.
3: examine
Simulation
resultsof fortheclusters
performancethe tilesize.
We
performancereconstruction
ofthe(A) Averaged
kernelwidth(# of voxels)
tile/region
size (#novelof voxels)
Increasing
the correlation
betweenhas(B) Averaged
little ofin5 clusters
of different
sizes
(from
5-35 voxels
per cluster).
0
0.1
0.3
3
5
7
9
11 13 15
5 10 15 20 25 30 355 40
algorithm
asafunctionoflocalsmoothingkernelwidth
andperformance
offunction
all voxels
algorithms.
onclusters
offluencedifferent
sizes
(from
per Basedcluster).
(A) Resultsof onthe theperformanceas a5-35
of increasing(A) Number of Clusters(B) Size of Clusters
0.2
Fig.
3:
Simulation
resultsoftheperformanceoftheplots,
it is We
clearexaminethat smoothChampagne
outperforms
thereconstructiontheofthe width
thesize.
local
smoothingkernelfrom
3 voxelsto121(A) Results
ofoftiletheperformance
as aperformance
functionoffor
increasing
0.1
allclusters
benchmarks.
Champagne
again
MxNE,
voxels
the
tile
size
was maintained
atoutperforms
6smoothingvoxels.
(B)
357
911 13 15
5 10 15 20 25 30 35 40
algorithm
asafunctionof
local
kernelwidt
5whileofdifferent
sizes
(from
5-35
voxelspercluster).
the
width
of
the
local
smoothingkernelfrom
3voxelsto121(A) Number of Clusters(B) Size of Clusterswhich
outperforms
MSP,
sLORETAanda Beamformer.
of theperformanceasperformancea functionof average
tile size of increasing(A)
Resultsof the
as
0.9
0.9
the
tile
size.
examineperformancefor sizereconstructi
voxelsResults
while
the
tile Wesize
wasmaintained
at 6function
voxels.
(B)
rangingfrom 4-132
voxels
per
tile.
The
results areaveraged
5 theInfluence
of
width
of
local
smoothingkernel
and
0.8
0.8
width
of
the
local
smoothingkernelfrom
3voxelsto
Resultsoftheperformanceasafunctionofaveragetile
size
5 clusters
sizes
(from
5-35
per121clu
over
50 of
simulations
atdifferent
each3data
andthe
the
error
bars
showvoxels
tiles of- Figure(A)point
shows
performanceof the0.7
0.7
voxels
while
the
tile
size
was
maintainedat6
voxels.
(B)
rangingfrom
4-132
voxelstile. The
resultsareaveragedthestandard
error.
novel
algorithm
asper
a performance
functionof localsmoothingkernel of incre
0.6
0.6
(A)
Results
of
the
as
a of
function
0.9
Resultsoftheperformanceasafunctionaverage
tile
size
0.5
0.5
over 50 simulations
each data
point andforthereconstruction
error bars show
width. Weatexamine
performanceof 5
0.8
the width
of
the
local
smoothingkernelfromvoxels t
rangingfrom
4-132
voxels
per
tile.
The
results
are3 averaged
0.4
0.4
clustersofdifferent
sizes
(from
5-35
voxelspercluster).
the standard error.
0.7
0.3
0.3
while
MCE
three
clusters
withseveralfocalerror
over
50localizes
simulations
atofeachdata
pointandkernel
theshowvoxelstheallwidth
tile
size
was
maintained
atbars6 voxelWe whileincreasethethe localsmoothing
from
0.2
0.2
0.6
point sources.
In tocontrast,
MSP lcalizes
all averaged
three clusters,
but wasthestandard
error.
3
voxels121
voxels
while
the
tile
size
Results of are
the performance
as a andfunctionof average til
0.1
0.1
reconstructions
groundadditional
0dB 3dB 6dB 9dB 0.5
12dB 15dB
0.1
0.3
0.5
0.7
0.9
maintained smootheraround 6 thanvoxels.
For truth
small cluster
sizes (of 5(C) SNR(D) Corr Between Clusters
thatare not 4-132
present
invoxels
theclusters
simulations
are produced.
ranging
from
per
tile.
The results
0.4
while sources
MCEvoxels)
localizes
all threewithseveralfocal are ave
performance
deteriorates
whensmoothingkernelBeamformer
sLORETA
MxNE
MSP
Champ
Smooth Champ Champagne can localize all three clusters but it estimates
0.3
point sources.
Inis
contrast,
MSP
all threeclusters,
buterror bars
overitsize
50
simulations
at lcalizes
each
data
point
and
the
increased.
However,
for threemoderateto large
clusters,
mostaccurately
localizes
clusters
whencompared
to00
activity that is more focal than the true spatial extent of the
0
0.2
Fig.
2:
Simulation
resultsofA
metricwith
four
different
reconstructionsaresmootherthangroundtruth
andadditional
performanceisinvarianttokernel
size
up
to
42.
Only
for
Fig.Fig.
2: 2:
Simulation
metric withwithfourfourdifferent
different
Simulationresults
results ofof A
A metric
thegroundtruth.
thestandard
error.
sources. Smooth Champagne outperforms the benchmarks as
0.1
configurations.
(A)Increasing
Increasingclusternumber;(B)Increasing
kernelforsize
ofclusters
121
voxelssources
thatveryarenot presentinthesimulations
areperformance
produced.
(A)
Increasingclusternumber;IncreasingdB 3dB configurations.
6dB configurations.
9dB 12dB 15dB
0.1
0.3
0.7 (B)(B)Increasing
0.9(A)
cluster0.5
number;3) large
Summary
simulations
it mostaccurately
localizes
threewhen does
comparedto(C) SNR(D) Corr
BetweenClusters intra-cluster
clustersize;(C)Increasing
SNR;(D)
Increasingstart
todeteriorate.
Figure
3
(B)
shows
the
performanceclustersize;(C)Increasing
Increasing SNR;
SNR; (D)Increasing
intra-cluster
Champagne
canlocalizeallthreeclusters
butitestimates
cluster
size;
(C)
(D)
Increasing
intra-cluster
the groundtruth.
As
can
be
seen
from
the
simulation
results
and analyses
correlation. The results are averaged over 50 simulations atofis thenovel
algorithm
as atruefunction
the tileofsize.
correlation.
The
are averaged
50
at at 3)amformer correlation.
sLORETA The
MxNEresults
MSP
Champ over
SmoothChamp
activitythat
morethan thespatialof extent
theSummary
for focal
simulations
50simulations
simulations
above,
SmoothChampagne
outperforms
alltheeach data pointresults
and the are
erroraveraged
bars show overthe standard
error.
We again
examine localizes
performancethree
for reconstructionof benchmark
5
it mostaccurately
clusters
whenpoint and0 the error bars show the standard error.
can
be
seen
from
the
simulation
results
andanalyses
Champagne
outperforms
thebenchmarks
as compa
eacheachdatadata
point
and theerror bars show the standard sources.
error. AsSmooth0sourcereconstruction
algorithms
withvariouscomplexdisclustersofdifferent
sizes
(from
5-35
voxelsperclus2: Simulationresults
results ofof A
fourfourdifferent
Simulation
A metric
metricwith
withdifferent
above,
Smooth localizes
Champagne
outperforms
allwhenthe benchmark
thegroundtruth.configurations.
it mostaccurately
threeclusters
compared
toevaluation
tributedsourceNext,
weextendtheter).
Weexamineperformanceasafunctionofaveragegurations.
(A)Increasing
Increasing cluster(B)Increasing
source reconstruction algorithms with various complex distions. (A)clusternumber;number;(B)Increasing3)
Summary
simulations
the groundoftruth.
theperformanceusing
realvoxels
MEGperandtileEEG
data.
tile
size
rangingfor
from
4-132
while
the
r size; (C) Increasing
SNR; (D) andIncreasing
intra-cluster
MxNE, sLORETA
Beamformer
show no change tributed source configurations. Next, we extend the evaluation
averagedsmoothingkernelsize
wasmaintained
aroundze;(C)
Increasing
SNR;(D)
Increasing
intra-cluster3)
Summary
simulations
As canforbe
seen from the simulation results and anation. The
averaged
50clusters
simulations
at MSP,
with areincreasingclustersize,
MxNE
outperforms
whileresults
MCE
localizes
all overthree
with several
focal
7
voxels.
Here,
we observethat across all cluster sizes,
on.
resultsover
50 simulations
atAsbut
Results
realdatasets
can beB.seen
fromof thesimulation
results and analyses
sLORETAand Beamformer.
above,
SmoothChampagne
outperforms
all the benchdataThepoint
andsources.
the are
error
bars
show
the lcalizes
standard
point
Inaveraged
contrast,
MSP
allerror.
three clusters,
performance deteriorates when the tile size increases. For
3
Effectsofincreasing
SNR
In
the
subplot(C)of
Figure2,
above,
SmoothChampagne
outperforms
allthebenchmarka pointreconstructions
and the errorarebars
showthanthegroundstandard
1)reconstruction
MEG- Somatosensory
EvokedField Paradigm
smoother
truth error.
and additionalsourcewithvarious
optimalperformance
thealgorithms
size ofthe smoothingkernel comple
versusis plotted. Allare
algorithms
source reconstruction
with various
dis- evoked
sources SNR
that arenotperformance
present in scorethe simulations
produced.
Figure 5 algorithms
shows the resultsof the complex
somatosensory
0.55

A PrimeA PrimeA Prime

0.5

3

7

0.75

0.65
0.6

0.7

0.55

cluster size

0.65

0.715

26

0.65

4

0.8

0.75

8

0.7

0.65
13

0.75

cluster size
0.7

23 5
48
(B) Averaged tile/region
size15(# of voxels)
0.6

0.65

25
35

0.55

0.8

5
15
25
35

0.5

121
542
(A) Averaged kernel
width (# 15
of voxels)
0.6

0.55

0.6

0.6

0.5

3

7

0.55

15

26

42

121

7

15

26

42

4

8

13

23

48

132

0.55 (B) Averaged tile/region size (# of voxels)(A) Averaged kernel width (# of voxels)

3

0.5

132

25
35

121

0.5

4

8

13

23

48

1

A PrimeA Prime

0.5

A Prime

0.8

0.75

0.8

A Prime

5
15
25
35

0.6

A Prime

0.65

0.9

0.9

0.85

cluster size

A Primecluster size 0.85

A Prime1

tributed source configurations. Next, we extend the evaltributed source
we extend theevaluation
Champagne can localize all three clusters but it estimates
field configurations.
response due toNext,
somatosensory
stimulipresented to aoftheperformance
using
real
MEG
and
EEG data.
real MEG
EEG data.
activity that is more focal than the true spatial extentofofthetheperformance
subject’s using
right index
finger,andaveragederived from a total ofsources. Smooth Champagne outperforms the benchmarks as 240 trials. A peak is typically seen ⇠50 ms after stimulation in72) EEG simulations: Figure 4 shows a representative
example of localization results for an EEG simulation. In thisconfiguration, the SNR, correlation of dipole activities within
the cluster and between clusters are 20 dB, 0.9 and 0.9, respectively. As is shown, the performance of all algorithms for EEG
simulation is almost the same as the performance for MEG
simulations. Compared with the ground truth, Beamformer
and sLORETA localize all three cluster activities but produce
blurred and spurious solutions, while MxNE localizes all threeclusters with several focal point sources. In contrast, MSP
lcalizes all three clusters, but reconstructions are smoother thanground truth and additional sources that are not present in thesimulations are produced. Champagne can localize all threeclusters but it estimates activity that is more focal than the true
spatial extent of the sources. Smooth Champagne outperforms
the benchmarks as it most accurately localizes three clusters
when compared to the ground truth.
3) Summary for simulations: As can be seen from thesimulation results and analyses above, Smooth Champagne
outperforms all the benchmark source reconstruction algorithms with various complex distributed source configurations.
Next, we extend the evaluation of the performance using real
MEG and EEG data.
B. Results of real datasets
1) MEG - Somatosensory Evoked Field Paradigm: Figure
5 shows the results of the somatosensory evoked field response
due to somatosensory stimuli presented to a subject’s rightindex finger, average derived from a total of 240 trials. A peakis 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. All algorithms canlocalize this activation to the correct area of somatosensory
cortex with focal reconstructions. MxNE and Champagne are
the most focal while sLORETA produces the most diffuse
reconstruction.

0.9

A Prime0.8

0.7

cluster size
5
15

0.6

25
35

0.5

3

7

15

26

42

121

(A) Averaged kernel width (# of voxels)
0.9

0.8

A Primeand the tile/region size should be smaller than the true
size of the clusters. If the kernel width exceeds the true
cluster size, then some deterioration in performance isobserved. If the tile/region size exceeds the true cluster
size, performance also deteriorates. However, in real data
and absence of ground truth, the true cluster size isunknown. Our analysis suggests that for cluster sizes
between 15-35 voxels, which are thought to the sizeof “blobs” typically reported in fMRI studies and may
represent realistic clustered brain activity, performanceis fairly uniform for a range of kernel sizes from 340 and for tile sizes between 4-25. This indicates somerobustness in algorithm performance. However, examining Figure 3 (A) and (B) indicates that robustness isconferred more by kernel size than by tile/region size,
especially for smaller clusters. In this context, if there isa mismatch between the kernel width and the tile/region
size, indeed performance will deteriorate due to the leadfield columns which are divergent such that the regional
lead-field cannot fit the scalp topography.

0.7

cluster size
5
15

0.6

25
35
0.5

4

8

13

23

48

132

(B) Averaged tile/region size (# of voxels)
Fig. 3: Simulation results of the performance of the novel
algorithm as a function of local smoothing kernel width andthe tile size. We examine performance for reconstruction of
5 clusters of different sizes (from 5-35 voxels per cluster).
(A) Results of the performance as a function of increasing
the width of the local smoothing kernel from 3 voxels to 121
voxels while the averaged tile size was maintained at 6 voxels.
(B) Results of the performance as a function of average tile
size ranging from 4-132 voxels per tile while the averagedsmoothing kernel size was maintained around 7 voxels. The
results are averaged over 50 simulations at each data point andthe error bars show the standard error.
2) MEG - Auditory Evoked Fields: The localization resultsfor AEF data from four subjects are shown in Figure 6. The
power of at each voxel in a 50-75 ms window around M100
peak is plotted for every algorithm. MxNE, Champagne andSmooth Champagne are able to consistently localize bilateralauditory activity for all subjects (shown in the last two columns
in Figure 6). The activity is in Heschl’s gyrus, which is thelocation of primary auditory cortex. The localization results of
Smooth Champagne are more spatially distributed than MxNE
and Champagne. Beamformer and sLORETA can find the two
auditory cortices only in one subject, whereas for the restof the subjects the activations are mostly biased towards thecenter of the head; suggesting that the correlation of bilateralauditory cortical activity greatly and negatively impacts theperformance of Beamformer and sLORETA. MSP can localize

8Cluster 1

Ground Truth

Beamformer

Right Cerebrum
Frontal Lobe (R)Middle Frontal
Gyrus (R)White Matter (R)Cluster 2

Right Cerebrum
Frontal Lobe (R)Middle Frontal
Gyrus (R)White Matter (R)Cluster 3

Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)sLORETA

MxNE

MSP

Right Cerebrum
Frontal Lobe (R)Middle Frontal
Gyrus (R)White Matter (R)Right Cerebrum
Frontal Lobe (R)Middle Frontal
Gyrus (R)White Matter (R)Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Champagne

Right Cerebrum
Frontal Lobe (R)Superior Frontal
Gyrus (R)Brodmann area
8 (R)Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Smooth Champagne

Right Cerebrum
Frontal Lobe (R)Sub-Gyral (R)White Matter(R)Right Cerebrum
Frontal Lobe (R)Middle Frontal
Gyrus (R)White Matter (R)Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Left Cerebrum
Frontal Lobe (L)Superior Frontal
Gyrus (L)Brodmann area
6 (L)Right Cerebrum
Lingual Gyrus (R)White Matter (R)Right Cerebrum
Occipital Lobe (R)Middle Occipital
Gyrus (R)White Matter (R)Fig. 4: A single EEG simulation example of the localization results with three clusters for six algorithms: Beamformer,
sLORETA, MxNE, MSP, Champagne and Smooth Champagne. The ground truth location of clusters are shown for comparison.
In this configuration, the SNR, correlation of dipole activities within the cluster and between clusters are 20 dB, 0.9 and 0.9,
respectively.
sLORETA

MxNE

Beamformer

sLORETA

MxNE

MSP

Champagne

Smooth Champagne

0.3

1000

0.2
0.2

0

Left Cerebrum
Parietal Lobe (L)Postcentral
Gyrus (L)White Matter (L)0.1
0Left Cerebrum
Parietal Lobe (L)Postcentral
Gyrus (L)White Matter (L)

500

0

-0.1
-0.2

-500

-0.2
-0.4

-0.3

-1000

Champagne

Smooth Champagne

Subject 4

MSP

Left Cerebrum
Parietal Lobe (L)Postcentral
Gyrus (L)Brodmann area
40 Matter (L)Subject 3

0.4

Subject 2

Subject 1

Beamformer

8

4

0

5000
Left Cerebrum
Parietal Lobe (L)Postcentral
Gyrus (L)Brodmann area
40 Matter (L)

0

-4

-8

1000

Left Cerebrum
Parietal Lobe (L)Postcentral
Gyrus (L)White Matter (L)5000Left Cerebrum
Parietal Lobe (L)Inferior Parietal
Lobule (L)

-500
-5000

-1000

Fig. 5: Sensory Evoked Field localization results. All six
algorithms localize to somatosensory cortical areas. Here we
set the threshold as half of the maximum value.

bilateral auditory activity but with some location bias anddiffuse activation.
3) MEG - Audio-Visual Evoked Fields: Figure 7 shows
results of the audio-visual evoked fields for Smooth Champagne. In subplot (A) and (B) we show the brain activations
associated with the auditory stimulus. Smooth Champagne isable to localize bilateral auditory activity in Heschl’s gyrus
in the window around the M100 peak, shown in the first row
of Figure 7. 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 Figure 7. SmoothChampagne is able to localize a source in the medial, occipitalgyrus with a peak around 150 ms. We plot the power inthe window around this peak and the time course of theFig. 6: Auditory evoked field (AEF) results with four subjects
for six algorithms: Beamformer, sLORETA, MxNE, MSP,
Champagne and Smooth Champagne. The results from bothChampagne and Smooth Champagne are shown in the last twocolumns, which outperform the other benchmark algorithms
shown in the first to fourth columns.

source marked with the cross hairs. Our algorithm can localize
a later visual response with a time course that has power
extending past 150 ms. Compared to the results obtained withChampagne [14], Smooth Champagne shows more diffuse
activation and smoother time courses.
4) EEG - Face-Processing task: In Figure. 8, we present
the results from using five algorithms on the face-processing
task EEG data set. Figure 8 shows the average power, M100
peak power, and M170 peak power at different rows separately.
MxNE is unable to localize brain activity and localizes one
focal activity for all conditions. Champagne is able to localize
brain activity with sparse peaks at visual areas and fusiform
gyrus [25], [26]. Smooth Champagne demonstrates spatiallydistributed activity in the appropriate corresponding locations.
Performance of benchmarks algorithms on this dataset can be9Intensity

B
Intensity

A

Intensity

D
Intensity

C

Fig. 7: Audio-Visual data localization results from SmoothChampagne. Smooth Champagne is able to localize a bilateralauditory response at 100 ms after the simultaneous presentation of tones and a visual stimulus. For bilateral auditory
activity, the results of locations and time courses are shownin (A), (B). Smooth Champagne can localize an early visual
response at 150 ms after the simultaneous presentation of tones
and visual stimulus shown in (C) and (D).
sLORETA

MxNE

Champagne

Smooth Champ

M 170

M 100

0-400 ms

Beamformer

Fig. 8: Face processing (EEG) localization results for fivealgorithms: Beamformer, sLORETA, MxNE, Champagne, andSmooth Champagne. The first row is the average power
mapping from 0 ms to 400 ms, the second and third rows
are for peak power activity at 100 ms and 170 ms separately.
found in [14].
V. D ISCUSSION
This paper derives a novel and robust Empirical Bayesianalgorithm for accurate reconstruction of distributed source activity from MEG and EEG data. The algorithm, which we referto as Smooth Champagne, readily handles a variety of configurations of distributed brain sources under high noise and interference conditions, especially correlated distributed sources
with complex spatial extent - a situation that commonlyarises even with simple imaging tasks. Smooth Champagne
is based on a principled cost function which maximizes aconvex lower bound on the marginal likelihood of the data andresulting update rules are fast and convergent. The algorithm
displays significant theoretical and empirical advantages overmany existing benchmark algorithms for electromagnetic brain
imaging. We compared Smooth Champagne’s performancewith benchmark algorithms such as Beamformer, sLORETA,
MxNE, Champagne and MSP on both simulated and real
data and demonstrated its superiority. Simulations of MEG
and EEG data explored algorithm’s performance for complex
source configurations with highly correlated time-courses andhigh levels of noise and interference. Smooth Champagne
outperforms all benchmarks, including Champagne, for distributed source clusters as it demonstrates more accurate
spatial extent reconstructions regardless of SNR or source
correlations. In general it is difficult to evaluate localization
algorithm performance with real data since the ground truthis unknown. For this reason, we chose real data sets that havewell-established patterns of brain activity. Performance onthese real MEG and EEG data sets demonstrates the superior
ability of Smooth Champagne to localize real brain distributed
brain activity across a variety of tasks (SEF, AEF, AV, andface-Processing) when compared to benchmarks.
Smooth Champagne builds upon our prior framework ofChampagne for sparse brain source reconstructions with two
key ideas: local spatial kernel smoothing and hyperparameter
tiling both of which contribute to algorithm performance.
Kernel size and tile/region size appear to independent affectperformance. Our analysis suggests that for optimal performance the size of the smoothing kernel and the tile/region
size should match the true size of the clusters. If the kernel
width exceeds the true cluster size, then some deterioration in
performance is observed. If the tile/region size exceeds the true
cluster size, performance also deteriorates. However, in real
data the true cluster sizes are unknown. Across simulations,
we find that with a setting for the kernel and tile size around
10 voxels, i.e. the number of regions is 10% of the number ofvoxels in the lead field, we are able to accurately reconstructdistributed sources that are even larger in size. In principle,
the kernel width parameter of the smoothing kernel and thetile/region size of the hyperparameter tiling can be estimatedfrom data by computing the model evidence as a function ofkernel width and the tile/region size, we relegate this effort tofuture work.
Related approaches for estimating distributed sources haveused Type-1 penalized likelihood methods. Many publications
have examined mixed-norm (MxNE) optimization for MEG,
typically L1 priors to impose sparsity in space and L2 priors
are used to impose smoothness in time [2], [3], [27], [28],[29], [30], [31]. For example, Haufe et al. [27] proposed Focal
Vector Field Reconstruction (FVR) by defining the regulating
penalty function, which renders the solutions unique, as a
global l1-norm of local l2-norms. In contrast, Zhu et al. [28]proposed variation and wavelet based sparse source imaging(VW-SSI) to better estimate cortical source locations andextents, which utilizes the L1-norm regularization method withthe enforcement of transform sparseness in both variation
and wavelet domains. Extensions have also been proposed
for time-frequency sparsity decompositions [3]. The latentparameters of both of these approaches are often estimatedby a Second-Order Cone Program (SOCP). More recently,
iterative reweighting strategies (IRES) have been used topenalize locations that are less likely to contain any source andcan provide reasonable information regarding the location andextent of the underlying sources [3], [32]. In contrast to these
Type-I likelihood estimation methods, in this paper we use
a Type-II likelihood estimation method. We have previously
shown that Type-I likelihood methods which make use ofnon-factorial, lead-field and noise-dependent priors can be10shown to have a dual with Type-II likelihood methods withcomparable performance [33]. However, in general, as wehave shown here, performance of Type-II likelihood estimation
methods yields superior results.
In the future, the Smooth Champagne framework can beincorporated into related hierarchical algorithms that we haveproposed recently, such as tree Champagne [16], that combines voxel and region level inferences for mixed source
configurations which may include both distributed and sparsebrain sources. The region-based variance model in SmoothChampagne is different from other approaches in several ways[34], [35]. First, in contrast to Bayesian model averaging(BMA), which implements multiple variance models and thenaverages the individual primary current density (PCD) for eachmodel using gradient descent methods, Smooth Champagne
uses a single generative model for all sources and uses empirical Bayesian inference to learn the regions distribution using
closed-form update rules and built-in model selection features[4]. Second, unlike greedy algorithms like those proposed in
Babadi et al. and Friston et al. [34], [36], which are approximation to Bayesian inference and are highly sensitive toinitialization and have the possibility of sub-optimal solutions,our closed form solutions with principled cost functions thathave convergence guarantees for Bayesian inference with nosubspace selection or pruning. Finally, Smooth Champagne
also has faster update rules based on convex-bounds on true
marginal likelihood of the data [4].
Future directions of our work are motivated by several
further extensions of the current algorithm. First, all priors in
Smooth Champagne are assumed to be Gaussian. However, thealgorithm could easily be modified to incorporate additional
sparse priors both across space, and especially across time.
Even within a Gaussian scale mixture modeling framework it
is possible to incorporate many sparse priors implicitly [37],[38] and we plan to explore algorithm performance for such
non-Gaussian extensions. Further more the current algorithm
makes no assumptions about temporal smoothness for sourcereconstructions. It is possible to include various forms ofpriors about temporal information within our framework ofcovariance component estimation. One such idea is to use temporal basis functions [39], and another is to use autoregressivemodels for source distributions over time. Such autoregressivemodels can also include spatiotemporal correlations both insource time-courses and in the time-course of the backgroundnoise and interference in the data [40], [41], [42]. Anotheridea that we plan to pursue is to incorporate forward modelingerrors by probabilistic modeling of the lead-field kernel within
our framework, an approach similar to related to the following
recent work [13], [43], [44].

A PPENDIX A

D ERIVATION OF THE COST FUNCTION
Here, we derive the expression for the cost function
F(Ω, Φ) shown in Eq. (10). Taking the logarithm of themarginal likelihood of the data p(Y|Ω), then we make theuse of the form [20] (pp.244)
p(Y, Z|Ω)log p(Y|Ω) = Ep(Z|Y) log
p(Z|Y)= Ep(Z|Y) [log p(Y|Z)] + Ep(Z|Y) [log p(Z|Ω)](A.1)−Ep(Z|Y) [log p(Z|Y)] .

Substitution of Eqs. (4), (5) and (6) into Eq. (A.1) results in
T
T

1 P
y(t) − GZ(t) Σ −1y(t) − GZ(t)
ε
T t=1
T
R
P
P
1
(A.2)
−
zn (t)T (Ωr I)−1 z(t)n
T r=1,n∈ζr t=1− log |Σy |,log p(Y|Ω) = −where Σy is the model data covariance matrix expressed asΣy = GΥ GT + Σ ε .

(A.3)

Since log |Σ y | is a concave function with respect to Υ [4],we can find an auxiliary parameter, Φr (r = 1, . . . , R), thatsatisfies the relationship
−

R
P
r=1

tr (Φr IΥ r ) + Φ0 ≤ − log |Σ y | ,(A.4)

where, Φ0 is a scalar term [45] [46]. Then we define a costfunction F(Ω, Φ), such that
T
T

1 P
y(t) − GZ(t)y(t) − GZ(t) Σ −1
ε
T t=1
R
T
P
P
1
zn (t)T (Ωr I)−1 zn (t)−(A.5)
T r=1,n∈ζr t=1
RP−
tr(Φr Ωr I) + Φ0 .

F(Ω, Φ) = −r=1

Eq. (A.4) guarantees that the relationshipF(Ω, Φ) ≤ log p(Y|Ω)),(A.6)

always holds. That is, the cost function F(Ω, Φ) forms a convex lower bound on the logarithm of the marginal likelihoodp(Y|Ω)). When we maximize F(Ω, Φ) with respect to Ω andΦ, such Ω also maximizes the marginal likelihood.
ACKNOWLEDGMENT
The authors would like to thank Susanne Honma, Danielle
Mizuiri, and Anne Findlay for collecting much of the MEG
data, and all members and collaborators of the Biomagnetic
Imaging Laboratory for their support. This work was supported
by NIH grants R01EB022717, R01DC013979, R01NS100440,
R01DC017696, and UCOP-MRP-17-454755.
R EFERENCES[1] D. P. Wipf, J. P. Owen, H. T. Attias, K. Sekihara, and S. S. Nagarajan,“Robust Bayesian estimation of the location, orientation, and time courseof multiple correlated neural sources using MEG,” NeuroImage, vol. 49,
pp. 641–655, 2010.
[2] D. Wipf and S. Nagarajan, “Iterative Reweighted l1 and l2 Methods for
Finding Sparse Solutions,” IEEE Journal of Selected Topics in Signal
Processing, vol. 4, no. 2, pp. 317–329, 2010.

11[3] A. Gramfort, D. Strohmeier, J. Haueisen, M. S. Hämäläinen, andM. Kowalski, “Time-frequency mixed-norm estimates: Sparse M/EEG
imaging with non-stationary source activations,” NeuroImage, vol. 70,
pp. 410–422, 2013.
[4] D. Wipf and S. Nagarajan, “A unified Bayesian framework for
MEG/EEG source imaging,” NeuroImage, vol. 44, no. 3, pp. 947–966,
2009.
[5] M. S. Hämäläinen and R. J. Ilmoniemi, “Interpreting measured magnetic fields of the brain: Estimates of current distributions,” Helsinki
University of Technology, Tech. Rep. TKK-F-A559, 1984.
[6] M. S. Hämäläinen and R. J. Ilmoniemi, “Interpreting magnetic fields ofthe brain: minimum norm estimates,” Medical & biological engineering& computing, vol. 32, no. 1, pp. 35–42, 1994.
[7] A. M. Dale, A. K. Liu, B. R. Fischl, R. L. Buckner, J. W. Belliveau,
J. D. Lewine, and E. Halgren, “Dynamic statistical parametric mapping:
Combining fMRI and MEG for high-resolution imaging of corticalactivity,” Neuron, vol. 26, pp. 55–67, 2000.
[8] R. D. Pascual-Marqui, “Standardized low resolution brain electromagnetic tomography (sLORETA): technical details,” Methods and Findings
in Experimental and Clinical Pharmacology, vol. 24, pp. 5–12, 2002.
[9] S. E. Robinson and D. F. Rose, “Current source image estimation by
spatially filtered MEG,” in Biomagnetism Clinical Aspects, M. Hokeet al., Eds. Elsevier Science Publishers, 1992, pp. 761–765.
[10] M. E. Spencer, R. M. Leahy, J. C. Mosher, and P. S. Lewis, “Adaptive
filters for monitoring localized brain activity from surface potential timeseries,” in Conference Record for 26th Annual Asilomer Conference on
Signals, Systems, and Computers, November 1992, pp. 156–161.
[11] B. D. Van Veen, W. Van Drongelen, M. Yuchtman, and A. Suzuki,“Localization of brain electrical activity via linearly constrained minimum variance spatial filtering,” IEEE Trans. Biomed. Eng., vol. 44, pp.
867–880, 1997.
[12] K. Sekihara and B. Scholz, “Generalized Wiener estimation of threedimensional current distribution from biomagnetic measurements,” in
Biomag 96: Proceedings of the Tenth International Conference on
Biomagnetism, C. J. Aine et al., Eds. New York: Springer-Verlag,
1996, pp. 338–341.
[13] S. A. H. Hosseini, A. Sohrabpour, M. Akçakaya, and B. He, “Electromagnetic Brain Source Imaging by Means of a Robust Minimum
Variance Beamformer,” IEEE Transactions on Biomedical Engineering,
vol. 65, no. 10, pp. 2365–2374, 2018.
[14] J. P. Owen, D. P. Wipf, H. T. Attias, K. Sekihara, and S. S. Nagarajan,“Performance evaluation of the champagne source reconstruction algorithm on simulated and real M/EEG data,” Neuroimage, vol. 60, no. 1,
pp. 305–323, 2012.
[15] C. Stahlhut, H. T. Attias, D. Wipf, L. K. Hansen, and S. S. Nagarajan,“Sparse spatio-temporal inference of electromagnetic brain sources,”in International Workshop on Machine Learning in Medical Imaging.
Springer, 2010, pp. 157–164.
[16] C. Cai, K. Sekihara, and S. S. Nagarajan, “Hierarchical multiscale
bayesian algorithm for robust MEG/EEG source reconstruction,” NeuroImage, vol. 183, pp. 698–715, 2018.
[17] H. Hallez, B. Vanrumste, R. Grech, J. Muscat, W. De Clercq, A. Vergult,
Y. D’Asseler, K. P. Camilleri, S. G. Fabri, S. Van Huffel et al., “Review
on solving the forward problem in EEG source analysis,” Journal ofneuroengineering and rehabilitation, vol. 4, no. 1, p. 46, 2007.
[18] S. S. Nagarajan, H. T. Attias, K. E. Hild, and K. Sekihara, “A probabilistic algorithm for robust interference suppression in bioelectromagnetic
sensor data,” Statistics in medicine, vol. 26, no. 21, pp. 3886–3910,
2007.
[19] K. Friston, L. Harrison, J. Daunizeau, S. Kiebel, C. Phillips, N. TrujilloBarreto, R. Henson, G. Flandin, and J. Mattout, “Multiple sparse priors
for the M/EEG inverse problem,” NeuroImage, vol. 39, no. 3, pp. 1104–
1120, 2008.
[20] K. Sekihara and S. S. Nagarajan, Electromagnetic brain imaging: A
Bayesian perspective. Berlin, Heidelber: Springer-Verlag, 2015.
[21] S. S. Dalal, J. Zumer, V. Agrawal, K. Hild, K. Sekihara, and S. Nagarajan, “NUTMEG: a neuromagnetic source reconstruction toolbox,”
Neurology & clinical neurophysiology: NCN, vol. 2004, p. 52, 2004.
[22] F. Darvas, D. Pantazis, E. Kucukaltun-Yildirim, and R. Leahy, “Mapping
human brain function with MEG and EEG: methods and validation,”
NeuroImage, vol. 23, pp. S289–S299, 2004.
[23] K. Sekihara, “Computing Resolution for Neuromagnetic Imaging
Systems,” J Comput Eng Inf Technol 5: 3. doi: http://dx. doi.
org/10.4172/2324, vol. 9307, p. 2, 2016.
[24] J. G. Snodgrass and J. Corwin, “Pragmatics of measuring recognition
memory: applications to dementia and amnesia.” Journal of Experimental Psychology: General, vol. 117, no. 1, p. 34, 1988.

[25] R. N. Henson, G. Flandin, K. J. Friston, and J. Mattout, “A Parametric
Empirical Bayesian framework for fMRI-constrained MEG/EEG sourcereconstruction,” Human brain mapping, vol. 31, no. 10, pp. 1512–1531,
2010.
[26] N. Kanwisher, J. McDermott, and M. M. Chun, “The fusiform face area:a module in human extrastriate cortex specialized for face perception,”
Journal of neuroscience, vol. 17, no. 11, pp. 4302–4311, 1997.
[27] L. Ding and B. He, “Sparse source imaging in electroencephalography
with accurate field modeling,” Human brain mapping, vol. 29, no. 9,
pp. 1053–1067, 2008.
[28] W. Ou, A. Nummenmaa, J. Ahveninen, J. W. Belliveau, M. S.
Hämäläinen, and P. Golland, “Multimodal functional imaging using
fMRI-informed regional EEG/MEG source estimation,” Neuroimage,
vol. 52, no. 1, pp. 97–108, 2010.
[29] M. Zhu, W. Zhang, D. L. Dickens, and L. Ding, “Reconstructing
spatially extended brain sources via enforcing multiple transform sparseness,” NeuroImage, vol. 86, pp. 280–293, 2014.
[30] S. Haufe, V. V. Nikulin, A. Ziehe, K.-R. Müller, and G. Nolte,“Combining sparsity and rotational invariance in EEG/MEG sourcereconstruction,” NeuroImage, vol. 42, no. 2, pp. 726–738, 2008.
[31] B. Cassidy and V. Solo, “Spatially Sparse, Temporally Smooth MEG
Via Vector l0,” IEEE transactions on medical imaging, vol. 34, no. 6,
pp. 1282–1293, 2015.
[32] A. Sohrabpour, Y. Lu, G. Worrell, and B. He, “Imaging brain source
extent from EEG/MEG by means of an iteratively reweighted edge
sparsity minimization (IRES) strategy,” NeuroImage, vol. 142, pp. 27–
42, 2016.
[33] D. P. Wipf, B. D. Rao, and S. Nagarajan, “Latent variable bayesianmodels for promoting sparsity,” IEEE Transactions on Information
Theory, vol. 57, no. 9, pp. 6236–6255, 2011.
[34] K. Friston, C. Chu, J. Mourão-Miranda, O. Hulme, G. Rees, W. Penny,
and J. Ashburner, “Bayesian decoding of brain images,” Neuroimage,
vol. 39, no. 1, pp. 181–205, 2008.
[35] N. J. Trujillo-Barreto, E. Aubert-Vázquez, and P. A. Valdés-Sosa,“Bayesian model averaging in EEG/MEG imaging,” NeuroImage,
vol. 21, no. 4, pp. 1300–1319, 2004.
[36] B. Babadi, G. Obregon-Henao, C. Lamus, M. S. Hämäläinen, E. N.
Brown, and P. L. Purdon, “A subspace pursuit-based iterative greedy hierarchical solution to the neuromagnetic inverse problem,” NeuroImage,
vol. 87, pp. 427–443, 2014.
[37] D. Wipf, “Bayesian methods for finding sparse representations,” Ph.D.
dissertation, University of California, San Diego, 2006.
[38] J. Palmer, B. D. Rao, and D. P. Wipf, “Perspectives on sparse Bayesianlearning,” in Advances in neural information processing systems, 2004,
pp. 249–256.
[39] J. M. Zumer, H. T. Attias, K. Sekihara, and S. S. Nagarajan, “Probabilistic algorithms for MEG/EEG source reconstruction using temporal basis
functions learned from data,” NeuroImage, vol. 41, no. 3, pp. 924–940,
2008.
[40] K. Liu, Z. L. Yu, W. Wu, Z. Gu, Y. Li, and S. Nagarajan, “Bayesian electromagnetic spatio-temporal imaging of extended sources with Markov
Random Field and temporal basis expansion,” NeuroImage, vol. 139, pp.
385–404, 2016.
[41] J. Daunizeau, J. Mattout, D. Clonda, B. Goulard, H. Benali, and J.-M.
Lina, “Bayesian spatio-temporal approach for EEG source reconstruction: conciliating ECD and distributed models,” IEEE Transactions on
Biomedical Engineering, vol. 53, no. 3, pp. 503–516, 2006.
[42] M. Fukushima, O. Yamashita, T. R. Knösche, and M.-a. Sato, “MEG
source reconstruction based on identification of directed source interactions on whole-brain anatomical networks,” NeuroImage, vol. 105, pp.
408–427, 2015.
[43] C. Stahlhut, M. Morup, O. Winther, and L. K. Hansen, “SOFOMORE:
Combined EEG source and forward model reconstruction,” in Biomedical Imaging: From Nano to Macro, 2009. ISBI’09. IEEE International
Symposium on. IEEE, 2009, pp. 450–453.
[44] S. T. Hansen, S. Hauberg, and L. K. Hansen, “Data-driven forward
model inference for EEG brain imaging,” NeuroImage, vol. 139, pp.
249–258, 2016.
[45] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul, “Anintroduction to variational methods for graphical models,” Machinelearning, vol. 37, no. 2, pp. 183–233, 1999.
[46] S. Boyd and L. Vandenberghe, Convex optimization.
Cambridge
university press, 2004.