NeuroImage 258 (2022) 119369

Contents lists available at ScienceDirect

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

Empirical Bayesian localization of event-related time-frequency neuralactivity dynamics
Chang Cai a,b,∗, Leighton Hinkley b, Yijing Gao b, Ali Hashemi c,d,e, Stefan Haufe c,f,
Kensuke Sekihara g,h, Srikantan S. Nagarajan b,∗∗aNational Engineering Research Center for E-Learning, Central China Normal University, Wuhan, China
Department of Radiology and Biomedical Imaging, University of California, San Francisco, CA 94143-0628, United StatescBerlin Center for Advanced Neuroimaging, Charité Universitätsmedizin Berlin, Berlin, GermanydMachine Learning Group, Electrical Engineering and Computer Science Faculty, Technische Universität Berlin, GermanyeInstitut für Mathematik, Technische Universität Berlin, GermanyfBernstein Center for Computational Neuroscience Berlin, Berlin, Germanyg
Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45 Yushima, Bunkyo-ku, Tokyo 113–8519, Japanh
Signal Analysis Inc., Hachioji, Tokyo, Japanba r t i c l e

i n f oKeywords:
Dynamics of neural activity
Brain source power changes
Five-dimensional neuroimaging
Electromagnetic brain imaging
Bayesian inferencea b s t r a c t
Accurate reconstruction of the spatio-temporal dynamics of event-related cortical oscillations across human brain
regions is an important problem in functional brain imaging and human cognitive neuroscience with magnetoencephalography (MEG) and electroencephalography (EEG). The problem is challenging not only in terms of
localization of complex source conﬁgurations from sensor measurements with unknown noise and interference
but also for reconstruction of transient event-related time-frequency dynamics of cortical oscillations. We recently proposed a robust empirical Bayesian algorithm for simultaneous reconstruction of complex brain source
activity and noise covariance, in the context of evoked and resting-state data. In this paper, we expand upon thisempirical Bayesian framework for optimal reconstruction of event-related time-frequency dynamics of regional
cortical oscillations, referred to as time-frequency Champagne (TFC). This framework enables imaging of ﬁvedimensional (space, time, and frequency) event-related brain activity from M/EEG data, and can be viewed as atime-frequency optimized adaptive Bayesian beamformer. We evaluate TFC in both simulations and several real
datasets, with comparisons to benchmark standards - variants of time-frequency optimized adaptive beamformers(TFBF) as well as the sLORETA algorithm. In simulations, we demonstrate several advantages in estimating timefrequency cortical oscillatory dynamics compared to benchmarks. With real MEG data, we demonstrate across
many datasets that the proposed approach is robust to highly correlated brain activity and low SNR data, and isable to accurately reconstruct cortical dynamics with data from just a few epochs.

1. Introduction
Noninvasive functional brain imaging has made a tremendous impact in improving our understanding of the human brain. Functional
magnetic resonance imaging (fMRI) has been the predominant modality for imaging the functioning brain since the middle of 1990s, as it can
produce images with high spatial resolution. However, due to relianceon the blood oxygenation response, fMRI lacks the temporal resolution
required to image the dynamic and oscillatory spatiotemporal patterns
associated with activities in the brain. In contrast to fMRI, magnetoencephalography (MEG) and electroencephalography (EEG) are other popular non-invasive techniques which measure electromagnetic ﬁelds em-∗∗∗

anating from the brain with excellent temporal resolution. This levelof resolution has clinical utility, as high temporal and frequency resolution is essential for epileptic waveform activity localization (Guggisberg
et al., 2008; Velmurugan et al., 2019). It has also shown to be essential for cognitive neuroscience imaging studies using M/EEG that examine dynamic brain networks subserving perception, action and cognition (Borghesani et al., 2020; Herman et al., 2013; Hinkley et al., 2011;
Ranasinghe et al., 2020).
Since M/EEG sensors are located outside the head, examining sensor
data does not provide suﬃcient information about which brain structures generated the activity of interest. To extract more precise information about underlying brain activities, it is essential to apply elec-

Corresponding author at: National Engineering Research Center for E-Learning, Central China Normal University, Wuhan, China.
Corresponding author at: Department of Radiology and Biomedical Imaging, University of California, San Francisco, CA 94143-0628, United States.
E-mail addresses: ccai@ccnu.edu.cn (C. Cai), sri@ucsf.edu (S.S. Nagarajan).

https://doi.org/10.1016/j.neuroimage.2022.119369.
Received 22 January 2022; Received in revised form 21 April 2022; Accepted 9 June 2022
Available online 11 June 2022.
1053-8119/© 2022 Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/)

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

tromagnetic brain source localization algorithms which are solutions
to the ill-posed biomagnetic inverse problem. Many inverse algorithms
have been proposed for solving the bioelectromagnetic inverse problem, i.e., for estimating the parameters of neural source activity from
MEG and EEG sensor data, including linear (van Drongelen et al., 1996;
Gross et al., 2003; Hämäläinen and Ilmoniemi, 1994; Schmidt, 1986;
Stoica and Nehorai, 1990; Van Veen et al., 1997) and non-linear methods (Friston et al., 2008; Gramfort et al., 2013; Wipf et al., 2010a), distributed (Pascual-Marqui, 2002) vs sparse (Gramfort et al., 2013; Wipf
et al., 2010a).
In functional neuroimaging studies, two distinct types of neuralsource activity have been examined – phase-locked and non-phaselocked event-related brain activity. Phase-locked activity is usually obtained in evoked-response paradigms by averaging activity across repeated events. Many algorithms have been proposed in the literaturefor examining phase-locked activity (Gramfort et al., 2013; Hämäläinen and Ilmoniemi, 1994; Pascual-Marqui, 2002; Uutela et al., 1999;
Wipf et al., 2010a). A key feature contributing to accurate estimationof phase-locked activity is that the power of phase-locked activity following an event is greater than a pre-event baseline period, and thisactivity is dominated by low frequencies (1–8 Hz). In contrast, nonphase locked activity is reconstructed by extracting signal magnitude(and phase) in single trials. Importantly, such event-related power ﬂuctuations can manifest as either an increase or a decrease in power.
Furthermore, power increases and decreases will be frequency-speciﬁc.
There have been only a few papers that have described algorithms for
accurate reconstructions of time-window and frequency speciﬁc eventrelated non-phase locked brain activity (Dalal et al., 2008a; Gross et al.,
2001; Limpiti et al., 2008).
In our previous work, we proposed an adaptive optimised time frequency Beamformer (TFBF) which enables accurate reconstruction ofﬁve-dimensional event-related non-phase locked brain activity by using
segmented ﬁltered data to obtain the weight matrix and is able to estimate the change in power from a control (i.e., baseline) time window toan active time window (Dalal et al., 2008a). However, the method suffers from sensitivity to inaccurate sample data covariance matrices thatare used to obtain the weight matrix, especially under high correlated
brain sources or limited time series data availability and low signal-tonoise ratio (SNR), for example, in cognitive paradigms. Although modiﬁcations to the Beamformer algorithm have been developed to address
the issue of correlated source reconstruction, they have met with varying success and are further limited by requiring a-priori information
on source location (Dalal et al., 2006), being limited to scalar sources
rather than vector sources (Brookes et al., 2006), or by being limited totwo sources due to computation time constraints (Diwakar et al., 2010;2011).
Recently, for phase-locked brain activity, we proposed a family of robust empirical Bayesian algorithms, e.g., Champagne (Cai et al., 2020;
Hashemi et al., 2021b; Sekihara and Nagarajan, 2015; Wipf et al.,2010b) and its extensions (Hashemi et al., 2021a; 2021c; Hashemi and
Haufe, 2018), which are able to reconstruct complex brain activity andlearn the noise covariance simultaneously, and which has been successfully applied in M/EEG source imaging. In this paper, we expand thisempirical Bayesian framework for reconstruction of event-related nonphase locked time-frequency dynamics of regional cortical oscillations.
We refer to this extension as time-frequency Champagne (TFC). ThisTFC framework enables imaging of ﬁve-dimensional (space, time, andfrequency) event-related brain activity from M/EEG data. The key idea
here that allows for reconstructions of event-related activity in bothtask and control periods is the use of sparse Bayesian reconstruction
algorithms on omnibus data, which includes both task and control timewindows (Robinson and Vrba, 1999; Singh et al., 2003; Vrba and Robinson, 2001).
Here, we evaluate TFC both in simulations and on several real
datasets, with comparisons to benchmark standards - variants of timefrequency optimized adaptive beamformers (TFBF) as well as the

sLORETA algorithm. In simulations, we demonstrate several advantages
in estimating time-frequency cortical oscillatory dynamics compared tobenchmarks. With real MEG data, we demonstrate across many datasetsthat the proposed approach is robust to highly correlated brain activity
and low SNR data, and is able to accurately reconstruct cortical dynamics with data from just a few epochs.
2. Methods
This section ﬁrst brieﬂy introduces the generative model for neuraldynamics and the problem of imaging dynamics of frequency-speciﬁcneural activity. Following this, we introduce the time-frequency Champagne algorithm (TFC) algorithm for estimation of neural dynamics. Finally, we show the relationship of TFC to TFBF.
2.1. The generative model
The sensor level data measured on or near the scalp can be modeledas:
𝒀 = 𝑳𝑺 + 𝜺 =

𝑁
∑
𝑛=1

𝑳𝑛 𝑺 𝑛 + 𝜺,(1)where 𝒀 = [𝒚 (1), ⋯ , 𝒚 (𝑇 )] ∈ ℝ𝑀×𝑇 is the measured sensor data, 𝑀 is thenumber of channels measured, and 𝑇 is the number of time samples.
𝑳 = [𝑳1 , … , 𝑳𝑁 ] ∈ ℝ𝑀×𝑑𝑐 𝑁 is the leadﬁeld matrix obtained from theforward model, and 𝑑𝑐 𝑁 = 𝑑𝑐 × 𝑁, which we assume to have 𝑑𝑐 orienta⊤
𝑑𝑐 𝑁×𝑇 is the unknown brain
tions for each voxel. 𝑺 = [𝑺 ⊤, … , 𝑺⊤𝑁] ∈ ℝ1𝑑activity. 𝑺 𝑛 = [𝒔1𝑛 , ⋯ , 𝒔𝑛𝑐 ]⊤ ∈ ℝ𝑑𝑐 ×𝑇 is the 𝑛-th voxel intensity. 𝑁 denotes𝑑the number of voxels under consideration and 𝑳𝑛 = [𝒍1𝑛 , ⋯ , 𝒍𝑛𝑐 ] ∈ ℝ𝑀×𝑑𝑐
denotes the leadﬁeld matrix for the 𝑛-th voxel. The 𝑖-th column of 𝑳𝑛
represents the signal vector that would be observed at the scalp given aunit current source or dipole at the 𝑛-th voxel with a ﬁxed orientation in
the 𝑖-th direction. The voxel dimension 𝑑𝑐 is usually set to 3. In eventrelated data analysis, we always segment the sensor data into control
and task periods which are expressed as 𝒀 𝐶 and 𝒀 𝑇 .
2.2. Imaging dynamics of frequency-speciﬁc neural activity
In many neuroscience paradigms, it is necessary to compute the dynamics of event-related non-phase locked neural activity for speciﬁc
experimental or cognitive tasks, which could be obtained by contrasting activity in multiple activation windows, across diﬀerent frequency
bands, or by comparison to activity during one or multiple baseline windows. The resulting contrast spectrograms are time frequency representation of non-phase locked event-related brain source activity.
In previous work, we have proposed an adaptive spatial ﬁltering
algorithm optimized for robust source time frequency reconstruction,
which maximizes the ultimate potential of M/EEG for ﬁve-dimensional
imaging of brain activity in space, time, and frequency and renders it
applicable for widespread studies of human cortical dynamics duringcognition (Dalal et al., 2008b). We directly computed spectrograms of
the source time series during a task time-period of interest by contrasting it with the spectrogram of a control time period where the brain
sources of interest are expected to not be active.
To quantify this task-induced activity, ﬁrst, the sensor data 𝒀 ispassed through a ﬁlter bank and then divided into control and task periods, which are denoted as 𝒀 𝐶 (𝑓 ) and 𝒀 𝑇 (𝑓 ) for 𝑓 -th band. Task periods
are then segmented into overlapping task windows, which are expressed
as 𝒀 𝑇 (𝑘, 𝑓 ) (𝑘 = 1, ⋯ , 𝐾, we assume the number of segmented task time
windows is 𝐾, the lower case 𝑘 refers to the index of the overlapping task
time window). To estimate the change from a control time window to a𝑘-th task time window at 𝑛-th voxel, we ﬁrst estimate the source powerfrom each windows at 𝑛-th voxel, denoted as 𝑃𝐶 (𝑛, 𝑓 ) and 𝑃𝑇 (𝑛, 𝑘, 𝑓 ),then the following equation to calculate the F-ratio of power contrast in
decibels is expressed as:
𝑃 (𝑛,𝑘,𝑓 )𝑭 (𝑛, 𝑘, 𝑓 ) = 10 log10 𝑃𝑇 (𝑛,𝑓 ) ,
𝐶

2

(2)

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Since MEG and EEG measure signals data outside brain, we need toestimate the underlying brain source activity from sensor data in bothtask and control time windows using inverse reconstruction algorithms.
Naive use of the above F-ratio to quantify task-induced brain activity in conjunction with sparse reconstruction algorithms can howeveroften lead to divide-by-zero errors for the most interesting voxels of interest, i.e. the voxels which are activated in the task-window but haveno activity in the control period resulting a zero in the denominator.
Moreover, if there are strong correlated brain sources, then beamformer
reconstruction of power in both task- and control-windows will be inaccurate, rendering the F-ratio calculation inaccurate as well. Here, we
propose a solution to avoid both of these problems with an adaptationof Champagne, a sparse Bayesian algorithm.

bands; then the ﬁltered active sensor data is segmented into overlapping active windows. Subsequently, TFC estimates the model data covariances of control, active and concatenated data using the Champagne
algorithm, and ﬁnally calculates the neural dynamics using Eq. (5).
2.4. Relationship of TFC to TFBF
Previously, we have proposed a TFBF algorithm to estimate brain
sources and their neural dynamics as (Dalal et al., 2008b):
𝑭 𝑇 𝐹 𝐵𝐹 (𝑛, 𝑘, 𝑓 )In our prior work, we introduced the robust empirical Bayesian algorithm Champagne, which is able to reconstruct complex brain activity
and estimate the source activation power (Cai et al., 2020). However,instead of estimating the Champagne source reconstruction estimate forpower in task and control time-windows separately, and to directly compute the power ratio at each voxel, here we propose a modiﬁed approach
for Eq. (2). Since Champagne is a sparse reconstruction algorithm, naive
estimation of the denominator, i.e. source power in control period, will
result in zero activity for many voxels and consequently result in divideby-zero errors for the calculation of the F-ratio in Eq. (2).
To circumvent this problem and to improve the robustness of the
F-ratio calculation, we estimate the power in task and control periods and propose using a spatial weight matrix derived from the omnibus data that includes both task and control time-windows 𝒀 𝑂 (𝑘, 𝑓 ) =[𝒀 𝐶 (𝑓 ), 𝒀 𝑇 (𝑘, 𝑓 )]. We then estimate omnibus hyperparameters 𝜶 𝑂 andthe well-regularized omnibus data covariance 𝚺𝑂 (𝑘, 𝑓 ), and use these
inferred parameters to estimate the source power in task and control
windows. The idea of using omnibus data has been used in conjunction with adaptive beamforming for MEG by others (Robinson and Vrba,
1999; Singh et al., 2003; Vrba and Robinson, 2001) and ensures robustness and avoids the divide-by-zero errors.
This modiﬁed time-frequency Champagne algorithm with noise
learning (Cai et al., 2021) is shown in Appendix A. Speciﬁcally, we obtain the brain source power from 𝑘-th task window and the control window at 𝑛-th voxel in 𝑓 -th frequency band 𝑃̂𝑇 (𝑛, 𝑘, 𝑓 ) and 𝑃̂𝐶 (𝑛, 𝑘, 𝑓 ),(3)

2 (𝑛, 𝑘, 𝑓 )tr (𝑳⊤ 𝚺−1 (𝑘, 𝑓 )𝚺 (𝑓 )𝚺−1 (𝑘, 𝑓 )𝑳 ),𝑃̂𝐶 (𝑛, 𝑘, 𝑓 ) = 𝛼𝑂
𝐶
𝑛
𝑛 𝑂
𝑂

(4)

−1
−1
tr (𝑳⊤𝑛 𝚺𝑂 (𝑘, 𝑓 )𝚺𝑇 (𝑘, 𝑓 )𝚺𝑂 (𝑘, 𝑓 )𝑳𝑛 )−1−1
tr (𝑳⊤𝑛 𝚺𝑂 (𝑘, 𝑓 )𝚺𝐶 (𝑓 )𝚺𝑂 (𝑘, 𝑓 )𝑳𝑛 ),tr (𝑾 ⊤ (𝑛,𝑘,𝑓 )𝑹𝑇 (𝑘,𝑓 )𝑾 (𝑛,𝑘,𝑓 )),tr (𝑾 ⊤ (𝑛,𝑘,𝑓 )𝑹𝐶 (𝑓 )𝑾 (𝑛,𝑘,𝑓 ))

⊤ −1
−1 .
𝑾 (𝑛, 𝑘, 𝑓 ) = 𝑹−1
𝑂 (𝑘, 𝑓 )𝑳𝑛 (𝑳𝑛 𝑹 𝑂 (𝑘, 𝑓 )𝑳𝑛 )

(7)

(8)

In the TFBF model above, the output power change for each voxel isdetermined by the lead ﬁeld matrix 𝑳𝑛 , the omnibus data covariance
𝑹𝑂 (𝑘, 𝑓 ), and the sample data covariances from task 𝑹𝑇 (𝑘, 𝑓 ) and control windows 𝑹𝐶 (𝑓 ).
The ﬁdelity of power change estimation is dramatically aﬀected by
the accuracy of the sample data covariance, as well as a variety of otherfactors such as SNR, source correlations, and the number of time samples 𝑇 . Since the TFBF algorithm requires that the covariance matrices be invertible, without regularization or with incorrect regularization for covariance matrix inversion, algorithm performance can degrade. Furthermore, if sources are highly correlated, the performanceof the estimate can be signiﬁcantly compromised, a well-known problem for adaptive beamformers. Finally, if the number of sample time
points 𝑇 is small, then the sample data covariance may deviate sharply
from the true covariance, leading to systematic estimation errors. Thisis partly because when 𝑇 is small, then even sources produced by an
uncorrelated generative model may appear correlated. Consequently,
if a more robust version of the sample data covariance were obtained,
it would greatly improve the power change estimation. The TFC algorithm we propose here contains a robust way to estimate the power ofour brain activity from control and task time window. Theoretically,
TFC can be shown to be a particular instantiation of a Bayesian modelbased TFBF, which are adaptive beamformers that use model covariances instead of the sample covariance of the data (Al-Shoukairi andRao, 2019; Pote and Rao, 2020). In the limit of inﬁnite data, the model
covariance matrix estimated by TFC will equal the sample covariance,
i.e. 𝑹 = 𝑳𝜶𝑳⊤ + 𝚲 = 𝜮 𝑦 . Therefore, we could replace 𝑹𝑂 (𝑘, 𝑓 ), 𝑹𝑇 (𝑘, 𝑓 )and 𝑹𝐶 (𝑓 ) with the corresponding model covariance matrices 𝜮 𝑂 (𝑘, 𝑓 ),
𝜮 𝑇 (𝑘, 𝑓 ), and 𝜮 𝐶 (𝑘, 𝑓 ) respectively in Eqs (7) and (8). Consequently,
the F-ratio of power contrast for TFBF can be expressed as:

where, 𝛼𝑂 (𝑛, 𝑘, 𝑓 ) is the prior variance of the 𝑛-th voxel for the omnibus
data, 𝚺𝑇 (𝑘, 𝑓 ), 𝚺𝐶 (𝑓 ), and 𝚺𝑂 (𝑘, 𝑓 ) are well-regularized model covariances for the task, control and omnibus data respectively which are robust to limited numbers of data samples. Consequently, when 𝛼𝑂 (𝑛, 𝑘, 𝑓 )is not equal to zero, this F-ratio as described above corresponds to voxels that are active in either the task or the control periods, which isexpressed as:𝑭 𝑇 𝐹 𝐶 (𝑛, 𝑘, 𝑓 ) = 10 log10

10 log10

where, 𝑃̂𝑇𝐵𝐹 (𝑛, 𝑘, 𝑓 ) and 𝑃̂𝐶𝐵𝐹 (𝑛, 𝑘, 𝑓 ) are the brain source power from 𝑘th task window and the control window at 𝑛-th voxel in 𝑓 -th frequency
band using Beamformer algorithm, 𝑹𝑇 (𝑘, 𝑓 ) ∈ ℝ𝑀×𝑀 and 𝑹𝐶 (𝑓 ) ∈ℝ𝑀×𝑀 are sample task and control covariances computed for each resulting time-frequency window. 𝑾 (𝑛, 𝑘, 𝑓 ) ∈ ℝ𝑀×𝑑𝑐 is the weight matrix
obtained with the omnibus sample data covariance, which is deﬁned as
𝑹𝑂 (𝑘, 𝑓 ).

2.3. Time-frequency Champagne algorithm

2 (𝑛, 𝑘, 𝑓 )tr (𝑳⊤ 𝚺−1 (𝑘, 𝑓 )𝚺 (𝑘, 𝑓 )𝚺−1 (𝑘, 𝑓 )𝑳 ),𝑃̂𝑇 (𝑛, 𝑘, 𝑓 ) = 𝛼𝑂
𝑇
𝑛
𝑛 𝑂
𝑂

=

𝑭 𝑇 𝐹 𝐵𝐹 (𝑛, 𝑘, 𝑓 )

=
=

−1
−1
tr (𝑳⊤𝑛 𝚺𝑂 (𝑘,𝑓 )𝚺𝑇 (𝑘,𝑓 )𝚺𝑂 (𝑘,𝑓 )𝑳𝑛 )−1−1
tr (𝑳⊤𝑛 𝚺𝑂 (𝑘,𝑓 )𝚺𝐶 (𝑓 )𝚺𝑂 (𝑘,𝑓 )𝑳𝑛 )𝑭 𝑇 𝐹 𝐶 (𝑛, 𝑘, 𝑓 ).

(9)

Therefore, TFC can be viewed as a particular instantiation of a Bayesianmodel beamformer that uses model covariance derived from sparse
Bayesian learning instead of the sample data covariance.

(5)

In contrast, while 𝛼𝑂 (𝑛, 𝑘, 𝑓 ) equals to zero, there is no activity in either
period, then the F-ratio can be expressed as:

3. Performance evaluation of simulated and real data

𝑭 𝑇 𝐹 𝐶 (𝑛, 𝑘, 𝑓 ) = 0.

3.1. Benchmark algorithms(6)

The speciﬁc steps of the time-frequency Champagne algorithm (TFC) aredetailed in Table 1.
In summary, to estimate the ﬁve-dimensional dynamics of neuralactivity, TFC ﬁrst ﬁlters the input sensor data into diﬀerent frequency

We chose to test the TFC algorithm against four representative
source localization algorithms from the literature: three adaptive spatial ﬁltering methods, namely time frequency optimized beamformer(TFBF) (Sekihara et al., 2001), Bayesian PCA based TFBF (bPCA_TFBF)
3

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369(Engemann and Gramfort, 2015; Woolrich et al., 2011), and Variational Bayesian Factor Analysis based TFBF (VBFA_TFBF) (Engemann
and Gramfort, 2015; Nagarajan et al., 2007), more details are provided
in Appendix B. In addition, one non-adaptive spatial ﬁltering method,sLORETA (Pascual-Marqui, 2002), is included. For sLORETA, we use thedefault setting deﬁned in the NUTMEG software where the regularization parameter of sLORETA is set as 1e-3 of the maximum eigenvalue of
the Gram matrix (Dalal et al., 2004).

time courses at frequencies of 8 Hz, 21 Hz and 42 Hz were synthesized
and placed randomly in the brain source space. We also randomized thephase of each source and varied them between trials. The brain sources
time series were windowed such that they represented event-related synchronization (ERS) activity and were not active simultaneously. The ﬁrst
source was active from 50 ms to 300 ms, while the second one was active from 350 ms to 550 ms, and the last source had a baseline activity
from -750 ms to 0 ms, then went to zero activity from 0 ms to 600 ms,
following which it returned to baseline activity from 600 ms to 1000 ms
to simulate event-related desynchronization (ERD) (Dalal et al., 2008a).
Then, we increased the number of seeded point sources to evaluate
the algorithm performance as a function of the number of sources. The
locations for the sources were randomly chosen and we also set someminimum distance between sources (at least 15 mm) and a minimum
distance from the center of the head (at least 15 mm). The generated
time courses were partitioned into control and active periods. The control period (-300 ms to 0 ms, sampled at 1200 Hz) contained only noise
and interfering brain activity. The active period of the same length (0 ms
to 300 ms, sampled at 1200 Hz) consisted of the activity of interest andinterference present in the control period. The SNR was ﬁxed at 5 dB.
The number of seeded point sources was increased from two to ten inincrements of two for 50 times tests per each conﬁguration. The interand intra source correlation strengths were set to be 0.95 and 0.95, formore details refer to Owen et al. (2012).
Last, we evaluated algorithm performance as a function of SNR. Simulations were performed at SNRs ranging from 0 dB to 20 dB in steps of
5 dB. The estimation of F-ratio performance was evaluated for ﬁve randomly seeded point sources with the same conﬁguration as describedabove for 50 tests. The inter and intra source correlation strengths were

3.2. Simulations conﬁgurations
Dipolar sources with variable orientation were used to generate thesimulated data. Damped sinusoidal time courses were created as point
source time activity and then projected to the sensors using the leadﬁeld matrix generated by the forward model. The brain volume was segmented into 8 mm voxels and a three-orientation forward lead ﬁeld was
calculated using a single spherical-shell model Hallez et al. (2007) implemented in NUTMEG (Hinkley et al., 2020a). The parcellation at thisresolution yields around 5000 voxels. To evaluate the robustness of the
proposed algorithm to noise, we randomly chose segments of real brain
noise from actual resting-state sensor recordings of ten human subjects.
These recordings are presumed to contain only spontaneous brain activity and sensor noise. SNR and correlations between voxel time courses
were varied to examine algorithm performance. SNR and time course
correlation are deﬁned in a standard fashion (Cai et al., 2018; Owen
et al., 2012).
We ﬁrst examined the accuracy of TFC with simulations. Fifty trials
were generated, spanning -750 ms to 1000 ms per trial with a sampling
rate of 1200 Hz and an SNR of 0 dB. Three point sources with sine waves

Fig. 1. A single example to compare the performanceof F-ratio estimation across methods. Shown at the top
are the peaks and corresponding spectrograms of thethree simulated sources. In the rows below the reconstruction results using TFBF, sLORETA, bPCA_TFBF,VBFA_TFBF and TFC are shown. In each of those panels, the crosshairs mark the spatiotemporal peak of
the F-ratio for the reconstructed source, with the corresponding spectrogram shown below it. The timefrequency window for which the spatial activation
map is superimposed on the MRI is highlighted in the
spectrogram.

4

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 2. Localization results of aggregate performance asa function of the number of seeded point sources and thesignal-to-noise ratio. The inter- and intra- source correlation were set to 0.95 and 0.95. The results are averaged
across 50 simulations for each conﬁguration and the error bars represent the standard error.

Fig. 3. Auditory evoked ﬁeld (AEF) results for two
subjects obtained for all available trials. The results
using TFC are shown in the right-most column andthe results obtained using the benchmark algorithms
are shown in the other four columns.

set to be 0.95 and 0.95. Again, we set the minimum distance between
sources and from the center of the head both at least 15 mm.

Next, we evaluated the performance of all algorithms on data for the
picture naming task. For the visual picture naming test, an image of an
object is projected onto a screen (99 trials) and subjects are instructedto name the pictured object verbally; more details can be found in ourprevious paper (Hinkley et al., 2020b). Again, for the frequency-domain
methods, in Beta band (12–30 Hz), the active window was set as -50 to
900 ms and the control window was -200 to 0 ms; in Gamma band (65–90 Hz), the active window was set as -200 to 300 ms and the control
window was -150 to -50 ms (Hinkley et al., 2011). Both datasets were
processed with a 50 ms window step size.
To evaluate the robustness of TFC to sample size, we evaluate algorithm performance when the number of trials was limited to 10, 20, 40
and 60 for ﬁnger movement and picturing name. Based on our previous work (Hinkley et al., 2019), for ﬁnger movement task, the active
window was set as follows: Beta band 12–30 Hz with active window
set to be -50 ms to 150 ms relative to the button press, with -950 ms
to -750 ms as the baseline; Gamma band 65–90 Hz with active window
set to be 0 ms to 100 ms relative to the button press, with -900 ms to 800 ms as the baseline. For picture naming task, we set the time window
and frequencies as follows, Beta band 12–30 Hz with active window set
to be -50 ms to 150 ms relative to the button press, with -950 ms to-750 ms as the baseline; Gamma band 65–90 Hz with active window
set to be 0 ms to 100 ms relative to the button press, with -900 ms to-800 ms as the baseline. Each test was performed 30 times with the speciﬁc trials themselves chosen as a random subset of all available trials(Hinkley et al., 2011). Since there was no existing ground truth for thereal tasks, to test robustness and stability of the results, we compare theestimated results with F-ratio estimated with whole trials by TFC during
the evaluation. It is worth noting that, even though we apply the TFC
as reference ”ground-truth” here, results are similar when we use TFBF
as the reference, because as the number of trials increases performanceof TFBF approaches that of TFC.

3.3. 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 sampling rate. The leadﬁeld for eachsubject was calculated in NUTMEG (Hinkley et al., 2020a) using asingle-sphere head model and an 8 mm voxel grid. Each column
was normalized to have a norm of unity. Three real MEG data sets
were used to evaluate the application of the algorithms: 1. Auditory
Evoked Fields (AEF); 2. Finger movement data; 3. Picture naming data.
The ﬁrst two data sets have been reported in our prior publications,and details about the datasets can be found in Wipf et al. (2010a),
Dalal et al. (2008a), Owen et al. (2012) for ﬁnger movement and picturenaming.
We collected 120 trials of AEF data from two subjects; more details
related to the task can be found in our prior papers (Cai et al., 2020;2018). We calculate the change in power from a control time window
to the active time window in an extended Alpha band ranging from 4
to 12 Hz (Lehtelä et al., 1997). We set the active window to range from
0 ms to 200 ms, and the baseline window to range from -250 ms to 50 ms. We also collected 96 trials of ﬁnger movement data from anothersubject. We tested the change in power in the Beta and Gamma bands
from a control (i.e., baseline) time window to an active time window.
For the frequency-domain methods, in Beta band (12–30 Hz), the active
window was set to range from -250 to 350 ms and the control window
was deﬁned to range from -950 to -750 ms; in the Gamma band (65–90 Hz), the active window was set as -200 to 300 ms and the control
window was -900 to -800 ms (Hinkley et al., 2019).
5

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 4. Shown above are F-ratios obtained for left index ﬁnger movement data using diﬀerent algorithms.
The change in power from a control time window
to an active time window was superimposed on subject’s brain for the Beta (ﬁrst row) and Gamma bands(second row).

3.4. Data processing and quantifying algorithm performance

curve is large, then the hit rate is higher compared to the false positive rate. Hit rate (𝐻𝑅 ) was calculated by dividing the number of hit
point sources by the true number of point sources. The false rate (𝐹𝑅 )was deﬁned by dividing the number of potential false positive voxelsby the total number of inactive voxels for each simulation. The details of the 𝐴′ metric calculation can be perused in our previous paper (Cai et al., 2019). The 𝐴′ metric ranges from 0 to 1, with higher
numbers reﬂecting better performance (Cai et al., 2018; Darvas et al.,
2004; Owen et al., 2012; Sekihara, 2016). To calculate the mean andvariance of the 𝐴′ metric, 50 simulations were conducted for eachconﬁguration.
In real data sets, since we did not know the actual brain activity
power change, to test the robustness and stability of the results, we compared estimated results with the F-ratio estimated from all trial epochs(120 epochs for AEF data, 96 epochs for ﬁnger movement, 99 epochs for
picture naming) by TFC for each subject and algorithm. We compared
the spatial correlation between estimated F-ratio for a smaller numberof trials with the F-ratio estimated for all trial epochs and the results forincreasing trial numbers were averaged across 30 simulations for eachconﬁguration.

We run TFC in a lattice of time-frequency windows. The original data
were ﬁrst passed through a bank of 200th-order ﬁnite impulse response(FIR) bandpass ﬁlters and subsequently split into overlapping temporalwindows with a step size of 25 ms for all bands. In our ﬁlter design,
we chose to follow traditional MEG/EEG power band deﬁnitions as best
as possible. We deﬁned Theta-Alpha band as 4 − 12 Hz with 300 ms
windows, Beta band 12 − 30 Hz with 200 ms windows, low Gamma 30 −55 Hz with 150 ms windows and high Gamma band 65 − 90 Hz with 100
ms windows.
To quantify the performance of TFC and other benchmark algorithms in simulations, in which we knew the ground truth, we usedfree-response receiver operator characteristics (FROC) which shows
the probability for detection of a true source in an image versus the
expected value of the number of false positive detections per image(Cai et al., 2018; Darvas et al., 2004; Owen et al., 2012; Sekihara,2016). Based on the FROC, we compute an 𝐴′ metric (Cai et al., 2019;
Snodgrass and Corwin, 1988), which is an estimate of the area under the FROC curve for each simulation. If the area under the FROC
6

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 5. Single example showing F-ratios for left index
ﬁnger movement in Beta band obtained for diﬀerent
algorithms as a function of the number of trials.

4. Results

the left column of Fig. 2. Here, the inter- and intra-source correlationcoeﬃcients were ﬁxed to 0.95 and SNR was ﬁxed to 5 dB. All algorithms present the same trend, showing decreasing performance underincreasing source count. TFC consistently produced better results than
TFBF, bPCA_TFBF, VBFA_TFBF and sLORETA in high correlation settings. bPCA_TFBF and VBFA_TFBF show similar performance, which asslightly inferior to that of TFBF, but both outperformed sLORETA.
Performance as a function of SNR is plotted for all algorithms in theright column of Fig. 2. Reconstruction performance was evaluated for
ﬁve randomly seeded point sources with both inter- and intra-source
correlation coeﬃcients of 0.95. Again, TFC consistently shows better
performance than other benchmarks in high correlation settings.
In summary, we demonstrate from simulated data that, while alladaptive beamforming based benchmarks outperform the non-adaptive
method - sLORETA, TFC consistently outperforms all benchmarks including adaptive beamformering based methods. Since bPCA_TFBF and
VBFA_TFBF show similar performance on the simulated data, we only
present the results of VBFA_TFBF on real data in the following.

4.1. Simulation results
In Fig. 1, we show a single example to compare the accuracy of the
F-ratio estimation for three simulated point sources. The TFBF was able
to provide good results and produce accurate localization for all threesources, the intuition was that the presence of jitter in the source timecourse across trials allowed Beamformers to reconstruct source power
across trials even if they were highly correlated within a single trial.
However, results bled over to other frequency bands and the estimatedduration of the activition was longer (greater temporal extent) than the
ground truth activation, especially for sources 1 and 2. For sLORETA,
the spatial extent of the solution was broader and more blurry for all
three simulated sources, with peaks on the periphery of the deﬁned volume of interest in each case. The SNR for sLORETA was also clearly
weaker than that of the TFBF solutions. bPCA_TFBF and VBFA_TFBF
show similar results as TFBF, and correctly localized all three simulated
sources, however, showing some frequency bleed-over and a general
overestimation of the temporal duration of the activations of sources 1
and 2. Finally, the results of TFC were superior to those of other benchmarks. TFC was able to localize all three sources, and showed minimal
frequency spill-over. The temporal durations of the activation estimatedby TFC were more tightly reﬂecting those of the ground truth comparedto other solutions, and produced the highest SNR (color-bar in Fig. 1) inthe spectrograms.
Fig. 2 shows the aggregate brain activity power change estimation
performance achieved in the simulations. The results of all algorithms
as a function of the number of seeded point sources are presented in

4.2. Results of real datasets
4.2.1. Auditory evoked ﬁeld paradigm
Fig. 3 shows auditory evoked ﬁeld (AEF) localization results obtainedusing all trials from two subjects. The power change from the control
time window to the active time window at each voxel around the M100
peak is plotted for each algorithm. TFBF, sLORETA and VBFA_TFBF
failed to generate convincing source solutions, whereas TFC is able tolocalize the expected bilateral brain activation showing focal reconstructions for both subjects. Speciﬁcally, the TFC estimates localize to Hes7

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 6. Single example showing F-ratios for left index
ﬁnger movement in high Gamma band obtained for
diﬀerent algorithms as a function of the number oftrials.

chl’s gyrus in the temporal lobe, the location of the primary auditory cortex (A1). These datasets demonstrate that TFC is superior to other source
solutions when localizing power changes of highly correlated multiple
sources from real data, which is known to be the case for sources in theauditory cortex.

Evaluation of reconstruction algorithm performance against increasing number of trials for left index ﬁnger movement data are plotted inFig. 7. We compared the spatial correlation between the F-ratios as estimated using diﬀerent algorithms and the F-ratio estimated by TFC using
all available trials. While increasing the number of trials increased theperformance of all algorithms, TFC outperformed the other solutions for
all numbers of trials in both Beta and Gamma bands.

4.2.2. Finger movement data
Fig. 4 shows F-ratios obtained for left index ﬁnger movements using
diﬀerent source estimation algorithms. The power change from a control
time window to an active time window in the Beta (12-30Hz; Fig. 4, toprow) and Gamma bands (65-90Hz; Fig. 4, second row) is superimposedon subject’s MRI. As expected, the characteristic Beta band power (12–30 Hz) decrease in contralateral sensorimotor cortex was observed by
all algorithms 0-200ms following movement onset. For sLORETA, thissolution was blurry (i.e. spatially disperse) and biased towards the edgeof the brain in voxels outside of sensorimotor cortex. Additionally, TFBF,VBFA_TFBF and TFC localized a focal high Gamma (65–90 Hz) peak insensorimotor cortex. This activity was found to be more spatially focal
and temporally bound to the movement. For high Gamma, the sLORETA
solution was spatially dispersed and not centered on a functionally relevant region for this task (sensorimotor cortex).
Figs. 5 and 6 illustrate the performance of each algorithm using real
(left index ﬁnger movement) data in Beta and high Gamma bands forincreasing numbers of trials. TFC is the only algorithm able to localize
a signiﬁcant Beta band power decrease and a high Gamma peak with as
few as ten trials. In contrast, VBFA_TFBF was only able to localize thesesources with at least 20 trials, TFBF only localized the same sources
with 40 trials, and sLORETA required 60 trials to identify these sources.
Overall, TFC outperformed all benchmarks in this dataset with a limited
number of trials.

4.2.3. Picture naming task paradigm
Fig. 8 shows the reconstructed power change from the picture naming task obtained by each algorithm from all available trials. TFBF,VBFA_TFBF and TFC were able to localize the visual cortex for the Beta
band, similar to our previous ﬁndings (Hinkley et al., 2011), while thelocalization of sLORETA was weighted/biased outside of the cortex andclose to the edge of the brain. In the high Gamma band, both TFBF
and TFC accurately localized increased power to the primary visual cortices (Hinkley et al., 2011), while other benchmarks localized sources
far away (dorsal) from these expected locations. In addition, the localization of TFC was closer to visual cortex and centered around 200ms,
conﬁrming the superior spatial and temporal resolution of TFC observedin the simulations.
F-ratio results obtained by diﬀerent algorithms for increasing numbers of trials of the picture naming task in Beta and Gamma bands areshown in Figs. 9 and 10. Similar to the motor data (Figs. 5 and 6), TFC
was able to localize Beta power changes to visual cortex with as few as
ten trials in a fashion more spatially and temporally focal than othersource solutions. Although TFBF and sLORETA perform adequately inreconstructing task-related eﬀects with a few number of trials, activation extents are more broad. Similarly, in high Gamma band, TFC wasable to localize power change accurately using all numbers of trials.
8

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 7. Obtained F-ratios as a function of the number oftrials for left index ﬁnger movement data. Spatial correlation was calculated between the F-ratios estimatedby diﬀerent algorithms and the ground truth, which was
deﬁned as the F-ratio estimated from all available trials
using TFC. The results were averaged across 30 simulations for each conﬁguration. Error bars show the standard error of the mean.

Fig. 8. Shown above are F-ratio estimates obtained
for the picture naming task using TFBF, sLORETA,
VBFA_Beam and TFC. The functional maps are superimposed on each individual subject’s structural
MRI scan. In each of the panels, the crosshairs markthe spatiotemporal peak for the reconstructed source,with the corresponding spectrogram shown below
it. The functional map plotted on the MRI corresponds to the time-frequency window highlighted in
the spectrogram.

9

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 9. A single example showing F-ratios obtained foreach algorithm on the picture naming task data in Betaband as a function of the number of trials.

Source localization for the other solutions were not as accurate in space
and time as they were for TFC.
Evaluation results for increasing numbers of trials of the picture naming task are shown in Fig. 11. Similar to what was demonstrated forthe motor data, performance of all algorithms increased with increasing
numbers of trials, and TFC outperformed all other benchmarks in boththe Beta and high Gamma band. Even though sLORETA showed higher
spatial correlation than TFBF and VBFA_TFBF for these datasets when
the number of trials was less than 40, this gain in performance was lostfor sLORETA when the number of trials increased to 40.
In summary, consistent with our expectations from the simulateddata, we were able to demonstrate that TFC consistently outperformed
other benchmarks in several real datasets. Adaptive beamformer based
benchmarks outperformed the non-adaptive method sLORETA, especially for complex brain source power change estimation, which wasconsistent to our previous ﬁndings.

of brain activity in various tasks even with few trials. TFC displayed signiﬁcant theoretical and empirical advantages over the existing benchmarks TFBF, bPCA_TFBF, VBFA_TFBF and non adaptive sLORETA. We
also demonstrate theoretically that the TFC algorithm can be viewed
as a time-frequency optimized beamformer algorithm where model covariances learnt from sparse Bayesian learning are used instead of thesample data covariance.
Simulations were developed to explore the reconstruction performance of complex source conﬁgurations with highly correlated timecourses, multiple point sources, and high levels of noise and interference. These simulations demonstrated that TFC outperformed all benchmarks since it showed less frequency bleed-over of activity and the temporal extents of the activations were accurately captured by TFC. In addition, TFC produces highest SNR (colorbar in Fig. 1) in spectrograms
and it is more robust to complex brain source reconstruction especiallywith high correlated time-courses, multiple point sources, and high levels of noise and interference as it show higher 𝐴′ . In general, it is diﬃcult to evaluate localization algorithm performance with real data since
the ground truth is unknown. For this reason, we chose real data setsthat have well-established patterns of brain activity (AEF, ﬁnger movement data and picture naming task). We also demonstrated that TFC
is able to localize the expected bilateral brain activation with focal reconstructions for all subjects for AEF data analysis. TFC was also able
to ﬁnd a characteristic Beta band power decrease in contralateral sensorimotor cortex and a focal high Gamma peak in sensorimotor cortex
which conﬁrmed results obtained by our existing robust ﬁve dimension
beamformer (Dalal et al., 2008a). In a picture naming task analysis, thelocalization quality of TFC was superior to that of other benchmarks,
since the localized voxel was more centered on the time-frequency decomposition and the spread of activation was more restricted. TFC thus

5. Discussion
This paper introduces a novel alternative empirical Bayesian method
to localize the time-frequency dynamics of non-phase locked event related cortical activity. Since the proposed algorithm explicitly expands
the robust source estimation algorithm, Champagne, we refer to thisalgorithm as Time-Frequency Champagne (TFC). TFC is able to obtainthe spatiotemporal dynamics of cortical oscillations across brain sourcesand enables the ﬁve-dimensional imaging of space, time, and frequency
activity in the brain, thereby it is applicable for wide ranging sets ofstudies of human cortical dynamics during cognition. TFC is able tohandle correlated brain activity and resolve many distinct and functionally relevant brain areas, and estimates the minor power changes
10

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Fig. 10. A single example showing F-ratios obtainedfor each algorithm on the picture naming task data in
Gamma band as a function of the number of trials.

Fig. 11. Evaluation results observed for increasing numbers of trials on the picture naming task data in the Betaand high Gamma bands. Spatial correlation was calculated between the F-ratios estimated by diﬀerent algorithms and ground truth, which was deﬁned as the Fratio estimated by TFC using all avilable trials. The results were averaged across 30 simulations for each conﬁguration. Error bars show the standard error of themean.

not only improved the spatial but also the temporal resolution, yielding
the highest output SNR.
Since spontaneous brain activity has a very low signal-to-noise ratio
against background activity, many trials are typically required to reconstruct evoked ﬁelds. In real data analyses, we also explored the ability of
our novel TFC to achieve such reconstructions with few trials. TFC wasable to robustly localize brain activity power changes with few trials in
the ﬁnger movement and picture naming tasks, which can be considereda breakthrough in electromagnetic brain imaging, showing a potential
to reduce scanning times for children and other groups of subjects who
cannot tolerate long periods of data collection.
Prior eﬀorts in model data covariance estimation have used pure statistical models like probabilistic principal components analysis (PPCA)
or factor analysis (FA) models. In our prior work and as benchmarks

here, we have demonstrated that such pure statistical model based covariance estimates do not yield robust performance. In contrast, we canview TFC as a hybrid method using a combined biophysical and statistical model that yields robust model covariances which can be usedinstead of sample covariance estimates in adaptive beamformers.
The robust performance of TFC enables accurate and automated estimation of event related non-phase locked activity from a very small
number of trials. Therefore, TFC can signiﬁcantly advance translational
neuroscience studies using MEG imaging. In clinical populations, suchas those with neurodevelopmental or neurodegenerative disorders, collecting data with minimal head movements for long scan times can beparticularly challenging. In these populations, it may be possible to use
TFC with data from shorter scans with fewer number of trials with minimal loss in reconstruction performance.
11

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

a convex function using Jesnen’s inequality (see recent publication
Hashemi et al., 2021b on MM for solving this class of problems).
Finally, in contrast to ours and the above mentioned approaches, theapproach proposed by Grova et al. (Chowdhury et al., 2013) uses a
maximum entropy non-linear framework for source reconstruction
and this framework does not have a direct interpretation within a
Bayesian framework. Appropriate modiﬁcations to these alternative
sparse methods could result in similar performance and need to be
demonstrated in future work. Related work also includes Gianvittorio
et al. who developed a Bayesian method for localization of multiple
dipole in the frequency domain (Luria et al., 2019). Future studies
should compare their approach with the proposed TFC.

Algorithm 1 Time Frequency Champagne Algorithm.
Input: Sensor data 𝒀 𝐶 , 𝒀 𝑇 and leadﬁeld matrix 𝑳 = [𝑳1 , … , 𝑳𝑁 ]Output: Dynamics of neural activity 𝑭
1: function TFC(𝒀 𝐶 , 𝒀 𝑇 , 𝑳)
Filter 𝒀 𝐶 and 𝒀 𝑇 into several frequency bands 𝒀 𝐶 (𝑓 ) and2:
𝒀 𝑇 (𝑓 ), 𝑓 = 1, ⋯ , 𝐹
for 𝑓 = 1, ⋯ , 𝐹 do
3:
Segment 𝒀 𝑇 (𝑓 ) into several overlapping windows 𝒀 𝑇 (𝑘, 𝑓 ),4:
(𝑘 = 1, ⋯ , 𝐾)[𝜮 𝐶 (𝑓 ), ∼]=Champagne(𝒀 𝐶 (𝑓 ), 𝑳)5:
for 𝑘 = 1, ⋯ , 𝐾 do
6:[𝜮 𝑇 (𝑘, 𝑓 ), ∼]=Champagne(𝒀 𝑇 (𝑘, 𝑓 ), 𝑳)
7:
𝒀 𝑂 (𝑘, 𝑓 ) = [𝒀 𝐶 (𝑓 ), 𝒀 𝑇 (𝑘, 𝑓 )]8:
[𝜮 𝑂 (𝑘, 𝑓 ), 𝜶 𝑂 ]=Champagne(𝒀 𝑂 (𝑘, 𝑓 ), 𝑳)9:
for 𝑛 = 1, ⋯ , 𝑁 do
10:
if 𝛼𝑂 (𝑛, 𝑘, 𝑓 ) ≠ 0 then
11:
𝑭

12:

𝑇𝐹𝐶

(𝑛, 𝑘, 𝑓 ) = 10 log10

Acknowledgment
The authors would like to thank Danielle Mizuiri and Anne Findlay
for collecting much of the MEG data, and all members and collaboratorsof the Biomagnetic Imaging Laboratory for their support. CC acknowledges support from National Natural Science Foundation of China grant
62007013 and Hubei Provincial Natural Science Foundation of China
under Grant 2021CFB384. AH acknowledges scholarship support from
the Machine Learning/Intelligent Data Analysis research group at Technische Universität Berlin, the Berlin Mathematical School (BMS) and the
Berlin Mathematics Research Center MATH+. SSN acknowledges support from NIH grants R01EB022717, R01DC013979, R01NS10-0440,
R01DC176960, R01DC017091, R01AG062196, UCOP-MRP-17-454755,
T32EB001631, DOD CDMRP Grant W81XWH1810741, and an industry
research contract from Ricoh MEG Inc. SH acknowledges funding fromthe European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No.
758985).

−1−1
tr (𝑳⊤𝑛 𝚺𝑂 (𝑘,𝑓 )𝚺𝑇 (𝑘,𝑓 )𝚺𝑂 (𝑘,𝑓 )𝑳𝑛 )else
𝑭 𝑇 𝐹 𝐶 (𝑛, 𝑘, 𝑓 ) = 0
14:
end if
15:
16:
end forend for17:
end for18:
19:
return 𝑭
20: end function−1−1
tr (𝑳⊤𝑛 𝚺𝑂 (𝑘,𝑓 )𝚺𝐶 (𝑓 )𝚺𝑂 (𝑘,𝑓 )𝑳𝑛 )

13:

21:
22: function Champagne(𝒀 , 𝑳)

Initialize
𝜶 = diag(𝛼1 𝑰 𝑑𝑐 ×𝑑𝑐 , ⋯ , 𝛼𝑁 𝑰 𝑑𝑐 ×𝑑𝑐 ), 𝚲 = diag(𝜆1 , ⋯ , 𝜆𝑀 )randomly
repeat24:
𝜮 𝒀 = 𝑳𝜶𝑳⊤ + 𝚲
25:
26:
for 𝑛 = 1, ⋯ , 𝑁 do
𝑺̄ 𝑛 = 𝛼̂ 𝑛 𝑳⊤
𝚺−1 𝒀
27:
√ 𝑛 𝒀
23:

28:
29:
30:
31:

𝛼̂ 𝑛 =

Credit authorship contribution statement
Chang Cai: Conceptualization, Methodology, Software, Validation,
Formal analysis, Investigation, Data curation, Writing – original draft,
Writing – review & editing, Visualization, Funding acquisition. Leighton
Hinkley: Methodology, Data curation, Writing – review & editing.
Yijing Gao: Conceptualization. Ali Hashemi: Writing – review & editing. Stefan Haufe: Conceptualization, Methodology, Writing – review& editing. Kensuke Sekihara: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Writing – original draft, Writing –review & editing. Srikantan S. Nagarajan: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Data curation, Writing– original draft, Writing – review & editing, Visualization, Supervision,
Project administration, Funding acquisition.

⊤tr (𝑺̄ 𝑛 𝑺̄ 𝑛 )−1
𝑇 tr (𝑳⊤𝑛 𝚺𝒀 𝑳𝑛 )end forfor 𝑚 = 1√
, ⋯ , 𝑀 do
𝜆̂ 𝑚 =[(𝒀 −𝑳𝑺̄ )(𝒀 −𝑳𝑺̄ )⊤ ]𝑚𝑚
𝑇 [𝜮 −1]
𝑌 𝑚𝑚

end foruntil = log |𝜮 𝒀 | + tr (𝒀 ⊤ 𝜮 −1
𝒀 ) converges
34:𝒀
return 𝜮𝒀 , 𝜶35:
36: end function
32:

33:

Appendix A. Brief summary of the Champagne algorithm
The underlying framework for Bayesian inversion algorithms proposed by Friston et al. (2008), Sato et al. (2004),
Belardinelli et al. (2012) et al. is similar to ours. These are all algorithms attempting to parameterize the inverse model in terms of
empirical Bayesian estimation of coeﬃcients of covariance components.
In the case of Friston et al. (2008). these covariance components are
based on spatial bases of the adjacency matrices. In the case of Satoet al. (2004), the covariance components are based on fMRI spatialpriors. These approaches can be thought of as estimating multiple
sparse priors (MSP) as investigated by Belardinelli et al. (2012).
Although our framework is identical to the above publications, there
are diﬀerences in the exact algorithm for optimization. In contrast to
Friston, Sato and Belardinelli et al. who use variations of Restricted
Maximum Likelihood, Expectation-Maximization, Fisher Scoring and
Greedy Methods with slower and poorer convergence properties, weuse a convex-bounding framework with better and faster convergence(see Wipf and Nagarajan, 2009 for details), which also correspondto Majorization Minimization algorithm based on majorization by

This section brieﬂy introduces the Champagne algorithm withnoise learning, for a detailed derivation, we refer to our previous
paper (Cai et al., 2020). Champagne algorithm with noise learning is formulated based on an empirical Bayesian schema and the
input are only the sensor data 𝒀 and leadﬁled matrix 𝑳. Champagne assumes the generative model (Eq. (1)) with the following
prior distributions 𝒔(𝑡) ∼  (𝒔(𝑡)|𝟎, 𝜶) and 𝜺(𝑡) ∼  (0, 𝚲), where 𝚲 =diag(𝜆1 , 𝜆2 , ⋯ , 𝜆𝑀 ) and 𝜶 is a 𝑑𝑐 𝑁 × 𝑑𝑐 𝑁 block diagonal matrix expressed as 𝜶 = diag(𝛼1 𝑰 𝑑𝑐 ×𝑑𝑐 , ⋯ , 𝛼𝑁 𝑰 𝑑𝑐 ×𝑑𝑐 ), 𝛼𝑛 𝑰 𝑑𝑐 ×𝑑𝑐 is a prior variance𝑑𝑐 × 𝑑𝑐 matrix of 𝒔𝑛 and 𝑰 𝑑𝑐 ×𝑑𝑐 is a 𝑑𝑐 × 𝑑𝑐 identity matrix. To estimate
the model data covariance, we need to maximize the marginal likelihood
expressed below:
 = log |𝜮 𝑌 | + tr (𝒀 ⊤ 𝜮 −1
𝑌 𝒀 ),(A.1)

where 𝜮 𝑌 = 𝑳𝜶𝑳𝑇 + 𝚲 is the model data covariance. The unknown parameters are the prior hyperparameters 𝜶 and 𝚲. Champagne utilizes aconvex bounding (Hashemi et al., 2021b; Jordan et al., 1999) on the12

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

𝑝(𝝐) =  (𝟎, 𝜎 −1 𝑰 ).

logarithm of the marginal likelihood (model evidence) to estimate these
unknown hyperparameters, which results in the following fast and convergent update rules:−1
𝑺̄ 𝑛 = 𝛼̂ 𝑛 𝑳⊤
𝑛 𝚺𝒀 𝒀 ,√𝛼̂ 𝑛 =𝜆̂ 𝑚 =⊤

tr (𝑺̄ 𝑛 𝑺̄ 𝑛 )−1
𝑇 tr (𝑳⊤𝑛 𝚺𝒀 𝑳𝑛 )√The prior on 𝑨 is an automatic relevance determination (ARD) prior:

(A.2)

,

[(𝒀 −𝑳𝑺̄ )(𝒀 −𝑳𝑺̄ )⊤ ]𝑚𝑚.
𝑇 [𝜮 −1]
𝑌 𝑚𝑚

𝑝(𝑨) =  (𝟎, (𝜎Ψ)−1 ),(A.4)

𝑨 = 𝑹𝒀 𝑿 (𝑹𝑿𝑿 + Ψ)−1 ,𝑿 = 𝜎[𝑨⊤ 𝜎𝑨 + 𝑀(𝑹𝑿𝑿 + Ψ)−1 + 𝑰 ]−1 𝑨𝑌 ,1
𝜎 = 𝐾𝑀tr (𝑹𝒀 𝒀 − 𝑨𝑹𝒀 𝑿 ),1
Ψ = diag[ 𝑀𝑨⊤ 𝜎𝑨 + (𝑹𝑿𝑿 + Ψ)−1 ],(B.5)

where, 𝑹𝒀 𝒀 is the sample data covariance of 𝒀 , 𝑹𝒀 𝑿 and 𝑹𝒀 𝑿 are
the expected data covariance of 𝒀 𝑿 and 𝑿𝑿 Sekihara and Nagarajan (2015). After convergence, the covariance derived from the bPCA
model is given by:𝐶𝐴𝚺𝑏𝑃= 𝑨𝑨⊤ + 𝜎𝑰 𝑀 .
𝒀

(B.6)

We note that the model used for bPCA in Eq. (B.4) can be considered as a
special case of the model used in Factor Analysis. In bPCA, the error, 𝝐 isassumed to be a homoscedastic noise, 𝝐 ∼  (0, 𝜎𝑰 𝑀 ): the noise varianceis the same for all variables, here all sensors. In contrast, FA assumesa heteroscedastic noise, 𝝐 ∼  (0, 𝚺𝜖 ), where 𝚺𝜖 = diag(𝜎1 , ⋯ , 𝜎𝑀 ): thenoise variance diﬀers between sensors. In this case, the update rule ofΣ𝜖 is expressed as 𝚺𝜖 = 𝐾1 diag(𝑹𝒀 𝒀 − 𝑨𝑹𝒀 𝑿 ). The estimated model data
covariance as delivered by FA is given by
𝚺𝐹𝒀 𝐴 = 𝑨𝑨⊤ + 𝚺𝜖 .

(B.7)

This estimate of the model data covariance (Eqs. (B.6) and (B.7)) canbe used to compute the beamformer weights and the estimated source
power, and ﬁnally the F-ratio of the power contrast. In this paper, we
named Time_Frequency Beamformer with model covariance learnt bybPCA and VBFA as bPCA_TFBF and VBFA_TFBF.

where 𝚺𝐶 (𝑓 ) and 𝚺𝑇 (𝑘, 𝑓 ) are well-regularized model covariances ofcontrol and task windows estimated by running Champagne algorithmwith inputs 𝒀 𝐶 (𝑓 ), 𝑳 and 𝒀 𝑇 (𝑘, 𝑓 ).

References
Al-Shoukairi, M., Rao, B.D., 2019. Sparse signal recovery using mpdr estimation. In:
ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal
Processing (ICASSP). IEEE, pp. 5047–5051.
Belardinelli, P., Ortiz, E., Barnes, G., Noppeney, U., Preissl, H., 2012. Source reconstruction
accuracy of MEG and EEG bayesian inversion approaches. PloS one 7 (12), e51985.
Borghesani, V., Hinkley, L.B., Ranasinghe, K.G., Thompson, M.M., Shwe, W., Mizuiri, D.,
Lauricella, M., Europa, E., Honma, S., Miller, Z., et al., 2020. Taking the sublexicalroute: brain dynamics of reading in the semantic variant of primary progressive aphasia. Brain 143 (8), 2545–2560.
Brookes, M.J., Stevenson, C.M., Barnes, G.R., Hillebrand, A., Simpson, M.I., Francis, S.T.,
Morris, P.G., 2006. Beamformer reconstruction of correlated sources using a modiﬁed
source model. NeuroImage 34. 1454–65
Cai, C., Diwakar, M., Chen, D., Sekihara, K., Nagarajan, S.S., 2019. Robust Empirical
Bayesian Reconstruction of Distributed Sources for Electromagnetic Brain Imaging.
IEEE Transactions on Medical Imaging 39 (3), 567–577.
Cai, C., Hashemi, A., Diwakar, M., Haufe, S., Sekihara, K., Nagarajan, S.S., 2020. Robust
estimation of noise for electromagnetic brain imaging with the Champagne algorithm.
NeuroImage 225, 117411.
Cai, C., Hashemi, A., Diwakar, M., Haufe, S., Sekihara, K., Nagarajan, S.S., 2021. Robust
estimation of noise for electromagnetic brain imaging with the Champagne algorithm.
NeuroImage 225, 117411.
Cai, C., Sekihara, K., Nagarajan, S.S., 2018. Hierarchical multiscale Bayesian algorithm
for robust MEG/EEG source reconstruction. NeuroImage 183, 698–715.
Chowdhury, R.A., Lina, J.M., Kobayashi, E., Grova, C., 2013. MEG source localization of
spatially extended generators of epileptic activity: comparing entropic and hierarchical bayesian approaches. PloS one 8 (2), e55969.
Dalal, S.S., Guggisberg, A.G., Edwards, E., Sekihara, K., Findlay, A.M., Canolty, R.T.,
Berger, M.S., Knight, R.T., Barbaro, N.M., Kirsch, H.E., Nagarajan, S.S., 2008. Five-dimensional neuroimaging: Localization of the time-frequency dynamics of cortical activity. NeuroImage 40, 1686–1700.
Dalal, S.S., Guggisberg, A.G., Edwards, E., Sekihara, K., Findlay, A.M., Canolty, R.T.,
Berger, M.S., Knight, R.T., Barbaro, N.M., Kirsch, H.E., et al., 2008. Five-dimensional
neuroimaging: localization of the time–frequency dynamics of cortical activity. Neuroimage 40 (4), 1686–1700.

Appendix B. Bayesian Principal Components Analysis (bPCA) andFactor Analysis (FA) estimates of model data covariance
Bayesian PCA (bPCA) corresponds to a multivariate Gaussian model
where a random vector can be expressed as a weighted linear combination of random components superimposed by independent noise. The
covariance can be decomposed as the sum of a low rank matrix and ascaled identity matrix. With this statistical model, standard PCA is transformed into a latent variable model such as VBFA.
To derive the bPCA model, we use an extension of the standard bPCA
model (Engemann and Gramfort, 2015; Nagarajan et al., 2007; Woolrich
et al., 2011). We assume 𝐾 represents the number of components and𝒀 a sample sensor time series. The M-dimensional vector 𝒀 is then obtained from a 𝐾-dimensional random vector 𝑿 which is linearly transformed by 𝐾 latent factors forming a matrix 𝑨 of size 𝑀 × 𝐾, and arandom noise vector 𝝐. The generative model used in bPCA is then expressed as:(B.1)

where we assume that both 𝑿 and 𝝐 are mutually independent random
vectors drawn from spherical multivariate Gaussian distributions, of respective dimensionalities 𝐾 and 𝑀. Without loss of generality, the priordistributions of 𝑿 and 𝝐 are assumed as:𝑝 ( 𝑿 ) =  ( 𝟎, 𝑰 ) ,(B.4)

where 𝜎𝑰 and 𝜎Ψ are the prior precision of noise 𝝐 and weight matrix
𝑨. Using Variational Bayes (VB) inference to approximate the unknown
parameters 𝑨, 𝑿 , 𝜎 and Ψ by maximizing the variational free energy,
the update rules are as follows,(A.3)

In short, the Champagne algorithm with noise learning ﬁrst initializes
the voxel variances and noise covariance with random values, then updates the voxel variances and the noise covariance using Eqs. (A.2),(A.3) and (A.4) until the cost function Eq. (A.1) converges. Finally, theChampagne algorithm with noise learning outputs the model data covariance 𝜮 𝑌 , brain source activity 𝑺̄ and estimated prior distributions
̂variance 𝜶.
For robust estimation of source power for control and time windows
using the Champagne algorithm, we start by ﬁltering and segmenting
the sensor data into control time and 𝐾 task time windows, to calculate the change from control time window 𝒀 𝐶 (𝑓 ) to 𝑘-th task time window 𝒀 𝑇 (𝑘, 𝑓 ). To ensure that the magnitude of sources are comparable between the task and control periods and decreases the likelihoodof resolving false sources, we run Champagne algorithm with the input of omnibus data that includes both task and control time-windows,
denoted as 𝒀 𝑂 (𝑘, 𝑓 ) = [𝒀 𝐶 (𝑓 ), 𝒀 𝑇 (𝑘, 𝑓 )] and leadﬁeld matrix 𝑳 to obtain omnibus prior hyperparamters 𝜶 𝑂 (𝑘, 𝑓 ) and model data covariance𝚺𝑂 (𝐾, 𝑓 ). We then use the omnibus parameters to estimate the source
power in control and task periods at 𝑛-th voxel (Eqs. (3) and (4) in themain paper) expressed as[]𝑃̂𝐶 (𝑛, 𝑘, 𝑓 ) = tr 𝐸𝑝(𝑺 𝐶 |𝒀 𝐶 ) [(𝑺 𝐶 )𝑛 (𝑺 ⊤𝐶 )𝑛 ]
(A.5)
−1
−1
2
= 𝛼𝑂
(𝑛, 𝑘, 𝑓 )tr (𝑳⊤𝑛 𝚺𝑂 (𝑘, 𝑓 )𝚺𝐶 (𝑓 )𝚺𝑂 (𝑘, 𝑓 )𝑳𝑛 ),[]𝑃̂𝑇 (𝑛, 𝑘, 𝑓 ) = tr 𝐸𝑝(𝑺 𝑇 |𝒀 𝑇 ) [(𝑺 𝑇 )𝑛 (𝑺 ⊤𝑇 )𝑛 ]
(A.6)
−1
−1
2
= 𝛼𝑂
(𝑛, 𝑘, 𝑓 )tr (𝑳⊤𝑛 𝚺𝑂 (𝑘, 𝑓 )𝚺𝑇 (𝑘, 𝑓 )𝚺𝑂 (𝑘, 𝑓 )𝑳𝑛 ),𝒀 = 𝑨𝑿 + 𝝐,

(B.3)

(B.2)
13

C. Cai, L. Hinkley, Y. Gao et al.

NeuroImage 258 (2022) 119369

Dalal, S.S., Sekihara, K., Nagarajan, S.S., 2006. Modiﬁed beamformers for coherent source
region suppression. IEEE Trans. Biomed. Eng. 53. 1357–63
Dalal, S.S., Zumer, J., Agrawal, V., Hild, K., Sekihara, K., Nagarajan, S., 2004. NUTMEG: aneuromagnetic source reconstruction toolbox. Neurology & clinical neurophysiology:
NCN 2004, 52.
Darvas, F., Pantazis, D., Kucukaltun-Yildirim, E., Leahy, R., 2004. Mapping human brain
function with MEG and EEG: methods and validation. NeuroImage 23, S289–S299.
Diwakar, M., Huang, M.X., Srinivasan, R., Harrington, D.L., Robb, A., Angeles, A., Muzzatti, L., Pakdaman, R., Song, T., Theilmann, R.J., Lee, R.R., 2010. Dual-core beamformer for obtaining highly correlated neuronal networks in meg. NeuroImage 54,
253–263.
Diwakar, M., Tal, O., Liu, T.T., Harrington, D.L., Srinivasan, R., Muzzatti, L., Song, T.,
Theilmann, R.J., Lee, R.R., Huang, M.X., 2011. Accurate reconstruction of temporal correlation for neuronal sources using the enhanced dual-core MEG beamformer.
NeuroImage 56. 1918–28
van Drongelen, W., Yuchtman, M., van Veen, B.D., van Huﬀelen, A.C., 1996. A spatial
ﬁltering technique to detect and localize multiple sources in the brain. Brain Topography 9, 39–49.
Engemann, D.A., Gramfort, A., 2015. Automated model selection in covariance estimationand spatial whitening of MEG and EEG signals. NeuroImage 108, 328–342.
Friston, K., Harrison, L., Daunizeau, J., Kiebel, S., Phillips, C., Trujillo-Barreto, N., Henson, R., Flandin, G., Mattout, J., 2008. Multiple sparse priors for the M/EEG inverse
problem. NeuroImage 39 (3), 1104–1120. doi:10.1016/j.neuroimage.2007.09.048.
Gramfort, A., Strohmeier, D., Haueisen, J., Hämäläinen, M.S., Kowalski, M., 2013.
Time-frequency mixed-norm estimates: Sparse M/EEG imaging with non-stationary
source activations. NeuroImage 70, 410–422.
Gross, J., Kujala, J., Hämäläinen, M., Timmermann, L., Schnitzler, A., Salmelin, R., 2001.
Dynamic imaging of coherent sources: studying neural interactions in the humanbrain. Proceedings of the National Academy of Sciences 98 (2), 694–699.
Gross, J., Timmermann, L., Kujala, J., Salmelin, R., Schnitzler, A., 2003. Properties of MEG
tomographic maps obtained with spatial ﬁltering. NeuroImage 19, 1329–1336.
Guggisberg, A.G., Kirsch, H.E., Mantle, M.M., Barbaro, N.M., Nagarajan, S.S., 2008. Fastoscillations associated with interictal spikes localize the epileptogenic zone in patientswith partial epilepsy. Neuroimage 39 (2), 661–668.
Hallez, H., Vanrumste, B., Grech, R., Muscat, J., De Clercq, W., Vergult, A., D’Asseler, Y.,
Camilleri, K.P., Fabri, S.G., Van Huﬀel, S., et al., 2007. Review on solving the forward
problem in EEG source analysis. Journal of neuroengineering and rehabilitation 4 (1),46.
Hämäläinen, M.S., Ilmoniemi, R.J., 1994. Interpreting magnetic ﬁelds of the brain: minimum norm estimates. Medical & biological engineering & computing 32 (1), 35–42.
Hashemi, A., Cai, C., Gao, Y., Ghosh, S., Müller, K.-R., Nagarajan, S., Haufe, S., 2021. Joint
Learning of Full-structure Noise in Hierarchical Bayesian Regression Models. bioRxiv
doi:10.1101/2021.11.28.470264.
Hashemi, A., Cai, C., Kutyniok, G., Müller, K.-R., Nagarajan, S., Haufe, S., 2021. Uniﬁcation of sparse Bayesian learning algorithms for electromagnetic brain imaging withthe majorization minimization framework. NeuroImage 239, 118309.
Hashemi, A., Gao, Y., Cai, C., Ghosh, S., Müller, K.R., Nagarajan, S.S., Haufe, S., 2021. Efﬁcient hierarchical Bayesian inference for spatio-temporal regression models in neuroimaging. In: Thirty-Fifth Conference on Neural Information Processing Systems.
Hashemi, A., Haufe, S., 2018. Improving EEG Source Localization through Spatio-temporal
Sparse Bayesian Learning. Submitted to 26th European Signal Processing Conference(EUSIPCO).
Herman, A.B., Houde, J.F., Vinogradov, S., Nagarajan, S.S., 2013. Parsing the phonologicalloop: activation timing in the dorsal speech stream determines accuracy in speech
reproduction. Journal of Neuroscience 33 (13), 5439–5453.
Hinkley, L.B., Dale, C.L., Cai, C., Zumer, J., Dalal, S., Findlay, A., Sekihara, K., Nagarajan, S.S., 2020. NUTMEG: Open Source Software for M/EEG Source Reconstruction.
Frontiers in Neuroscience 14, 710.
Hinkley, L.B., Dale, C.L., Luks, T.L., Findlay, A.M., Bukshpun, P., Pojman, N., Thieu, T.,
Chung, W.K., Berman, J., Roberts, T.P., et al., 2019. Sensorimotor cortical oscillations
during movement preparation in 16p11. 2 deletion carriers. Journal of Neuroscience
39 (37), 7321–7331.
Hinkley, L.B., De Witte, E., Cahill-Thompson, M., Mizuiri, D., Garrett, C., Honma, S.,
Findlay, A., Gorno-Tempini, M.L., Tarapore, P., Kirsch, H.E., et al., 2020. Optimizing
Magnetoencephalographic Imaging Estimation of Language Lateralization for Simpler
Language Tasks. Frontiers in Human Neuroscience 14.
Hinkley, L.B., Nagarajan, S.S., Dalal, S.S., Guggisberg, A.G., Disbrow, E.A., 2011. Cortical
temporal dynamics of visually guided behavior. Cerebral Cortex 21 (3), 519–529.
Jordan, M.I., Ghahramani, Z., Jaakkola, T.S., Saul, L.K., 1999. An introduction to variational methods for graphical models. Machine learning 37 (2), 183–233.

Lehtelä, L., Salmelin, R., Hari, R., 1997. Evidence for reactive magnetic 10-hz rhythm inthe human auditory cortex. Neuroscience letters 222 (2), 111–114.
Limpiti, T., Van Veen, B.D., Attias, H.T., Nagarajan, S.S., 2008. A spatiotemporal framework for estimating trial-to-trial amplitude variation in event-related MEG/EEG. IEEE
Transactions on Biomedical Engineering 56 (3), 633–645.
Luria, G., Duran, D., Visani, E., Sommariva, S., Rotondi, F., Sebastiano, D.R., Panzica, F.,
Piana, M., Sorrentino, A., 2019. Bayesian multi-dipole modelling in the frequency
domain. Journal of neuroscience methods 312, 27–36.
Nagarajan, S.S., Attias, H.T., Hild, K.E., Sekihara, K., 2007. A probabilistic algorithm
for robust interference suppression in bioelectromagnetic sensor data. Statistics in
Medicine 26. 3886–910
Owen, J.P., Wipf, D.P., Attias, H.T., Sekihara, K., Nagarajan, S.S., 2012. Performance
evaluation of the champagne source reconstruction algorithm on simulated and real
M/EEG data. Neuroimage 60 (1), 305–323.
Pascual-Marqui, R.D., 2002. Standardized low resolution brain electromagnetic tomography (sLORETA): technical details. Methods and Findings in Experimental and ClinicalPharmacology 24, 5–12.
Pote, R.R., Rao, B.D., 2020. Robustness of sparse Bayesian Learning in Correlated Environments. In: ICASSP 2020-2020 IEEE International Conference on Acoustics, Speechand Signal Processing (ICASSP). IEEE, pp. 9100–9104.
Ranasinghe, K.G., Petersen, C., Kudo, K., Srivatsan, S., Beagle, A.J., Mizuiri, D., Findlay, A., Houde, J.F., Rankin, K., Rabinovici, G.D., et al., 2020. Alpha-frequency synchronization deﬁcits during life predict postmortem neuroﬁbrillary tangle burden inalzheimer’s disease: Biomarkers: Leveraging postmortem collections to validate neuroimaging. Alzheimer’s & Dementia 16, e045351.
Robinson, S., Vrba, J., 1999. Functional neuroimaging by synthetic aperture magnetometry (sam). Recent Advances in Biomagnetism,In: T. Yoshimoto, M. Kotani, S. Kuriki,
H. Karibe and N. Nakasato, Eds. 302–305.
Sato, M.-a., Yoshioka, T., Kajihara, S., Toyama, K., Goda, N., Doya, K., Kawato, M.,
2004. Hierarchical Bayesian estimation for MEG inverse problem. NeuroImage 23(3), 806–826.
Schmidt, R.O., 1986. Multiple emitter location and signal parameter estimation. IEEE
Trans. Antenn. Propagat. 34, 276–280.
Sekihara, K., 2016. Computing Resolution for Neuromagnetic Imaging Systems. J Comput
Eng Inf Technol 5: 3. doi: http://dx. doi. org/10.4172/2324 9307, 2.
Sekihara, K., Nagarajan, S.S., 2015. Electromagnetic brain imaging: A Bayesian perspective. Springer-Verlag, Berlin, Heidelber.
Sekihara, K., Nagarajan, S.S., Poeppel, D., Marantz, A., Miyashita, Y., 2001. Reconstructing spatio-temporal activities of neural sources using an MEG vector beamformer technique. IEEE Trans. Biomed. Eng. 48, 760–771.
Singh, K.D., Barnes, G.R., Hillebrand, A., 2003. Group imaging of task-related changes
in cortical synchronisation using nonparametric permutation testing. Neuroimage 19(4), 1589–1601.
Snodgrass, J.G., Corwin, J., 1988. Pragmatics of measuring recognition memory: applications to dementia and amnesia. Journal of Experimental Psychology: General 117 (1),34.
Stoica, P., Nehorai, A., 1990. MUSIC, maximum likelihood, and cramer-rao bound: further
results and comparisons. IEEE Trans. Acoust. Speech Signal Process. 38, 2140–2150.
Uutela, K., Hämäläinen, M., Somersalo, E., 1999. Visualization of magnetoencephalographic data using minimum current estimate. NeuroImage 10, 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.
Velmurugan, J., Nagarajan, S.S., Mariyappa, N., Mundlamuri, R.C., Raghavendra, K.,
Bharath, R.D., Saini, J., Arivazhagan, A., Rajeswaran, J., Mahadevan, A., et al., 2019.
Magnetoencephalography imaging of high frequency oscillations strengthens presurgical localization and outcome prediction. Brain 142 (11), 3514–3529.
Vrba, J., Robinson, S.E., 2001. Signal processing in magnetoencephalography. Methods
25 (2), 249–271.
Wipf, D., Nagarajan, S., 2009. A uniﬁed bayesian framework for MEG/EEG source imaging. NeuroImage 44 (3), 947–966.
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.
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 (1), 641–655.
Woolrich, M., Hunt, L., Groves, A., Barnes, G., 2011. MEG beamforming using Bayesian
PCA for adaptive data covariance matrix regularization. Neuroimage 57 (4),1466–1479.

14