Quantifying Time-Varying Sources in Magnetoencephalography -- A Discrete Approach

TL;DR

This paper introduces a novel discrete dynamic framework for localizing and quantifying time-varying brain sources in MEG data, addressing non-linear relationships and providing comprehensive temporal source distribution insights.

Contribution

The paper develops a new discrete dynamic approach for estimating source location distributions in MEG, overcoming limitations of existing linear methods and enabling full temporal source characterization.

Findings

The method accurately tracks source location changes over time.

It provides detailed source distribution insights from entire measurement sets.

The approach outperforms traditional methods in simulation and real data.

Abstract

We study the distribution of brain source from the most advanced brain imaging technique, Magnetoencephalography (MEG), which measures the magnetic fields outside the human head produced by the electrical activity inside the brain. Common time-varying source localization methods assume the source current with a time-varying structure and solve the MEG inverse problem by mainly estimating the source moment parameters. These methods use the fact that the magnetic fields linearly depend on the moment parameters of the source, and work well under the linear dynamic system. However, magnetic fields are known to be non-linearly related to the location parameters of the source. The existing work on estimating the time-varying unknown location parameters is limited. We are motivated to investigate the source distribution for the location parameters based on a dynamic framework, where the posterior distribution of the source is computed in a closed form discretely. The new framework allows us not only to directly approximate the posterior distribution of the source current, where sequential sampling methods may suffer from slow convergence due to the large volume of measurement, but also to quantify the source distribution at any time point from the entire set of measurements reflecting the distribution of the source, rather than using only the measurements up to the time point of interest. Both a dynamic procedure and a switch procedure are proposed for the new discrete approach, balancing estimation accuracy and computational efficiency when multiple sources are present. In both simulation and real data, we illustrate that the new method is able to provide comprehensive insight into the time evolution of the sources at different stages of the MEG and EEG experiment.

Tables

Table 1. Table 1: Forward procedure of the forward-backward algorithm.

Aim: Calculation of ${\alpha }_{tk}\left(𝚯\right)$, $1\le t\le T$, $1\le k\le K$.

Input: Parameter $𝚯$, and discretization ${\left\{V_k\right\}}_{k=1}^K$.

1. Compute $\text{P}\left(v_{1k}=1|𝚯\right)$, $\text{P}(v_{tk}=1|v_{t-1,l}=1,𝚯)$, $2\le t\le T$, and $\text{P}\left({𝐘}_t|v_{tk}=1,𝚯\right)$,

$1\le t\le T$.

2. Compute

$c_1\left(𝚯\right):=p\left({𝐘}_1|𝚯\right)\approx {\sum }_{k=1}^K\text{P}\left({𝐘}_1|v_{1k}=1,𝚯\right)\text{P}\left(v_{1k}=1|𝚯\right)$,

${\alpha }_{1k}\left(𝚯\right)=\text{P}\left({𝐘}_1|v_{1k}=1,𝚯\right)\text{P}\left(v_{1k}=1|𝚯\right)/c_1\left(𝚯\right)$, for $1\le k\le K$.

3. For $t=2,\dots ,T$, compute ${\alpha }_{tk}\left(𝚯\right)$ by using ${\alpha }_{t-1,l}\left(𝚯\right)$,

for $1\le k\le K$.

Output: ${\left\{{\alpha }_{tk}\left(𝚯\right)\right\}}_{t=1,k=1}^{T,K}$ and ${\left\{c_t\left(𝚯\right)\right\}}_{t=1}^T$

Table 2. Table 2: Backward procedure of the forward-backward algorithm.

Aim: Calculation of ${\beta }_{tk}\left(𝚯\right)$, $1\le t\le T$, $1\le k\le K$.
Input: Parameter $𝚯$, discretization ${\left\{V_k\right\}}_{k=1}^K$, and ${\left\{c_t\left(𝚯\right)\right\}}_{t=1}^T$ from the forward
procedure.
1. Compute $\text{P}\left(v_{1k}=1|𝚯\right)$, $\text{P}(v_{tk}=1|v_{t-1,l}=1,𝚯)$, $2\le t\le T$, and $\text{P}\left({𝐘}_t|v_{tk}=1,𝚯\right)$,
$1\le t\le T$.
2. Initialize ${\beta }_{Tk}\left(𝚯\right)=1$, for $1\le k\le K$.
3. For $t=T-1,\dots ,1$, compute ${\beta }_{tk}\left(𝚯\right)$ by using ${\beta }_{t+1,l}\left(𝚯\right)$ and $c_{t+1}\left(𝚯\right)$,
${\beta }_{tk}\left(𝚯\right)=1/c_{t+1}\left(𝚯\right){\sum }_{l=1}^K{\beta }_{t+1,l}\left(𝚯\right)\text{P}\left({𝐘}_{t+1}\right|v_{t+1,l}=1,𝚯)\text{P}(v_{t+1,l}=1|v_{tk}=1,𝚯)$,
for $1\le k\le K$.
Output: ${\left\{{\beta }_{tk}\left(𝚯\right)\right\}}_{t=1,k=1}^{T,K}$

Table 3. Table 3: Closed form for updates of estimated parameters at each iteration in the EM algorithm.

Update ${𝚯}^{\left(j\right)}=\{{𝝁}_0^{\left(j\right)},{𝚺}_0^{\left(j\right)},A^{\left(j\right)},{𝐛}^{\left(j\right)},{𝚺}^{\left(j\right)},{𝐕}^{\left(j\right)}\}$, given discretization ${\left\{V_k\right\}}_{k=1}^K$,
intermediate posterior distributions ${\left\{{\xi }_{tk}\left({𝚯}^{\left(j-1\right)}\right)\right\}}_{t=1,k=1}^{T,K}$ and ${\left\{{\eta }_{t-1,l}^{t,k}\left({𝚯}^{\left(j-1\right)}\right)\right\}}_{t=2,k,l=1}^{T,K}$,
and constant moment parameter $𝐪$.
${𝐝}_k:={({𝐜}_k^T,{𝐪}^T)}^T$,
${𝝁}_0^{\left(j\right)}={\sum }_{k=1}^K{\xi }_{1k}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_k$,
${𝚺}_0^{\left(j\right)}={\sum }_{k=1}^K{\xi }_{1k}\left({𝚯}^{\left(j-1\right)}\right)\left({𝐝}_k-{𝝁}_0^{\left(j\right)}\right){\left({𝐝}_k-{𝝁}_0^{\left(j\right)}\right)}^T$,
$A^{\left(j\right)}=[\left({\sum }_{t=2}^T{\sum }_{k=1}^K{\xi }_{tk}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_k\right)\left({\sum }_{t=2}^T{\sum }_{l=1}^K{\xi }_{t-1,l}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_l^T\right)-(T-1)\cdot$
$\left({\sum }_{t=2}^T{\sum }_{k=1}^K{\sum }_{l=1}^K{\eta }_{t-1,l}^{tk}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_k{𝐝}_l^T\right)]\cdot [\left({\sum }_{t=2}^T{\sum }_{l=1}^K{\xi }_{t-1,l}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_l\right)\cdot$
$\left({\sum }_{t=2}^T{\sum }_{l=1}^K{\xi }_{t-1,l}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_l^T\right)-(T-1)\left({\sum }_{t=2}^T{\sum }_{l=1}^K{\xi }_{t-1,l}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_l{𝐝}_l^T\right)]{}^{-1}$,
${𝐛}^{\left(j\right)}=\left({\sum }_{t=2}^T{\sum }_{k=1}^K{\xi }_{tk}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_k-A^{\left(j\right)}{\sum }_{t=2}^T{\sum }_{l=1}^K{\xi }_{t-1,l}\left({𝚯}^{\left(j-1\right)}\right){𝐝}_l\right)/\left(T-1\right)$,
${𝚺}^{\left(j\right)}=\frac{1}{T-1}\left({\sum }_{t=2}^T{\sum }_{k=1}^K{\sum }_{l=1}^K{\eta }_{t-1,l}^{tk}\left({𝚯}^{\left(j-1\right)}\right)\left({𝐝}_k-A^{\left(j\right)}{𝐝}_l-{𝐛}^{\left(j\right)}\right){\left({𝐝}_k-A^{\left(j\right)}{𝐝}_l-{𝐛}^{\left(j\right)}\right)}^T\right)$,
${𝐕}^{\left(j\right)}=\left({\sum }_{t=1}^T{\sum }_{k=1}^K{\xi }_{tk}\left({𝚯}^{\left(j-1\right)}\right)\left({𝐘}_t-𝐁\left({𝐝}_k\right)\right){\left({𝐘}_t-𝐁\left({𝐝}_k\right)\right)}^T\right)/T$.

Table 4. Table 4: EM algorithm.

Initialize ${𝚯}^{\left(0\right)}$. For $j=1,2,\dots$,

1. E-step. Calculate the posterior distribution ${\xi }_{tk}\left({𝚯}^{\left(j-1\right)}\right)$, for $1\le t\le T$, $1\le k\le K$,

and the joint posterior distribution ${\eta }_{t-1,l}^{tk}\left({𝚯}^{\left(j-1\right)}\right)$, for $2\le t\le T$, $1\le k,l\le K$.

2. M-step. Maximize the expected complete data log-likelihood defined in (2.16),

${𝚯}^{\left(j\right)}={argmax}_{𝚯}Q\left(𝚯|{𝚯}^{\left(j-1\right)}\right)$.

Table 5. Table 5: Source simulation in Case 1: The location parameter 𝐩 t subscript 𝐩 𝑡 {\bf p}_{t} of the source is expressed in terms of Cartesian coordinates ( x 𝑥 x (cm), y 𝑦 y (cm), z 𝑧 z (cm)) and is allowed to vary. The moments and strength parameter 𝐪 t subscript 𝐪 𝑡 {\bf q}_{t} is fixed during simulations.

${𝝁}_0={({𝐩}_0^T,{𝐪}_0^T)}^T$	${(-2,1,5,3,3,3)}^T$
${𝚺}_0=\text{diag}({\sigma }_{0,1}^2,{\sigma }_{0,2}^2,\dots ,{\sigma }_{0,6}^2)$	diag$(0.0225,0.0225,0.0225,{10}^{-4},{10}^{-4},{10}^{-4})$
$A=\text{diag}(a_1,a_2,\dots ,a_6)$	$\text{diag}(0.75,0.8,0.9,1,1,1)$
$𝐛={(b_1,b_2,\dots ,b_6)}^T$	${(0.75,-0.5,0.25,0,0,0)}^T$
$𝚺=\text{diag}({\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_6^2)$	diag $(0.25,0.25,0.25,{10}^{-4},{10}^{-4},{10}^{-4})$
$𝐕=\text{diag}({\sigma }^2,{\sigma }^2,\dots ,{\sigma }^2)$	diag$(6.25\ast {10}^{-5},6.25\ast {10}^{-5},\dots ,6.25\ast {10}^{-5})$
Number of time points	$100$

Table 6. Table 6: Comparison of mean errors for MLEs from the EM algorithm with different procedures. Standard deviations computed based on four repetitions are shown in parentheses.

	Case 1		Case 2
	dynamic EM	non-dynamic EM	dynamic EM	non-dynamic EM
			with switch procedure	with switch procedure
${\Vert \hat{A}-A\Vert }_{\text{max}}$	0.1293 (0.0078)	0.1642 (0.0128)	0.3346 (0.0577)	0.4074 (0.0578)
${\Vert \widehat{𝐛}-𝐛\Vert }_{\text{max}}$	0.2095 (0.0257)	0.3110 (0.1165)	1.7887 (0.1788)	2.3282 (0.4476)

Table 7. Table 7: Source simulation in Case 2: The location parameter 𝐩 t subscript 𝐩 𝑡 {\bf p}_{t} of the sources is expressed in terms of Cartesian coordinates ( x 𝑥 x (cm), y 𝑦 y (cm), z 𝑧 z (cm)) and is allowed to vary. The moments and strength parameter 𝐪 t subscript 𝐪 𝑡 {\bf q}_{t} is fixed during simulations.

${𝝁}_{0,1}={({𝐩}_{0,1}^T,{𝐪}_{0,1}^T)}^T$	${(1,1,5,3,3,3)}^T$
${𝝁}_{0,2}={({𝐩}_{0,2}^T,{𝐪}_{0,2}^T)}^T$	${(-1,2,4,3,3,3)}^T$
${𝚺}_{0,1}=\text{diag}({\sigma }_{0,1}^2,{\sigma }_{0,2}^2,\dots ,{\sigma }_{0,6}^2)$	$\text{diag}(0.01,0.01,0.01,{10}^{-4},{10}^{-4},{10}^{-4})$
${𝚺}_{0,2}=\text{diag}({\sigma }_{0,1}^2,{\sigma }_{0,2}^2,\dots ,{\sigma }_{0,6}^2)$	$\text{diag}(0.01,0.01,0.01,{10}^{-4},{10}^{-4},{10}^{-4})$
$A_1=\text{diag}(a_{1,1},a_{1,2},\dots ,a_{1,6})$	$\text{diag}(0.5,0.8,0.9,1,1,1)$
$A_2=\text{diag}(a_{2,1},a_{2,2},\dots ,a_{2,6})$	$\text{diag}(0.45,0.75,0.85,1,1,1)$
${𝐛}_1={(b_{1,1},b_{1,2},\dots ,b_{1,6})}^T$	${(2,-1,0.25,0,0,0)}^T$
${𝐛}_2={(b_{2,1},b_{2,2},\dots ,b_{2,6})}^T$	${(1.8,-0.8,0.5,0,0,0)}^T$
${𝚺}_1=\text{diag}({\sigma }_{1,1}^2,{\sigma }_{1,2}^2,\dots ,{\sigma }_{1,6}^2)$	diag$(0.25,0.25,0.09,{10}^{-4},{10}^{-4},{10}^{-4})$
${𝚺}_2=\text{diag}({\sigma }_{2,1}^2,{\sigma }_{2,2}^2,\dots ,{\sigma }_{2,6}^2)$	diag$(0.25,0.25,0.09,{10}^{-4},{10}^{-4},{10}^{-4})$
$𝐕=\text{diag}({\sigma }^2,{\sigma }^2,\dots ,{\sigma }^2)$	diag$(6.25\ast {10}^{-5},6.25\ast {10}^{-5},\dots ,6.25\ast {10}^{-5})$
Number of time points	$100$

Equations

$${\mathbf{J}}_n^{\text{p}}(r)$$

$${\mathbf{B}}_n(r)$$

$$Y_{l,t}$$

$$B_l({\mathbf{J}}_t^{\text{p}})$$

$$H_l({\mathbf{J}}_t^{\text{p}})$$

$${\mathbf{Y}}_t$$

$${\mathbf{J}}_t^{\text{p}}$$

$${\mathbf{J}}_{t,n}^{\text{p}}$$

$$\displaystyle p({\bf J}^{\text{\tiny p}}_{t}|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$})$$

$$\displaystyle p({\bf J}^{\text{\tiny p}}_{1}|\mbox{\boldmath$\Theta$})\approx\mbox{P}({\bf J}^{\text{\tiny p}}_{1}\in V_{k^{\ast}}|\mbox{\boldmath$\Theta$})$$

$$\displaystyle p({\bf J}^{\text{\tiny p}}_{t}|{\bf J}^{\text{\tiny p}}_{t-1},\mbox{\boldmath$\Theta$})\approx\Pi_{k=1}^{K}\Pi_{l=1}^{K}{\mbox{P}(v_{tk}=1|v_{t-1,l}=1,\mbox{\boldmath$\Theta$})}^{v_{tk}v_{t-1,l}},$$

$$\displaystyle p({\bf Y}_{t}|{\bf J}^{\text{\tiny p}}_{t},\mbox{\boldmath$\Theta$})\approx\Pi_{k=1}^{K}{\mbox{P}({\bf Y}_{t}|v_{tk}=1,\mbox{\boldmath$\Theta$})}^{v_{tk}},$$

$$\displaystyle\mbox{P}({\bf J}^{\text{\tiny p}}_{t}\in V_{k}|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$})$$

$$\displaystyle\hat{\mbox{\boldmath$\Theta$}}$$

$$\displaystyle\ell(\mbox{\boldmath$\Theta$},\mathcal{Y}_{T})$$

$$\displaystyle p(\mathcal{Y}_{T},\mathcal{J}^{\text{\tiny p}}_{T}|\mbox{\boldmath$\Theta$})$$

$$\displaystyle\ell(\mbox{\boldmath$\Theta$},\mathcal{Y}_{T})$$

$$=:$$

$$\displaystyle\ell(\mbox{\boldmath$\Theta$},\mathcal{Y}_{T})$$

$$\displaystyle\ell(\mbox{\boldmath$\Theta$},\mathcal{Y}_{T})=L(\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$}),\mbox{\boldmath$\Theta$},\mathcal{Y}_{T})=\log p(\mathcal{Y}_{T}|\mbox{\boldmath$\Theta$}).$$

$$\displaystyle\mbox{P}(v_{tk}=1|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$}^{(j-1)})=\text{E}(v_{tk}|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$}^{(j-1)})=\sum_{\mathcal{J}^{\text{\tiny p}}_{T}}\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$}^{(j-1)})v_{tk},$$

$$\displaystyle\mbox{P}(v_{t-1,l}=1,v_{tk}=1|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$}^{(j-1)})$$

$$\displaystyle\sum_{\mathcal{J}^{\text{\tiny p}}_{T}}\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath$\Theta$}^{(j-1)})\log p(\mathcal{Y}_{T},\mathcal{J}^{\text{\tiny p}}_{T}|\mbox{\boldmath$\Theta$})=:Q(\mbox{\boldmath$\Theta$}|\mbox{\boldmath$\Theta$}^{(j-1)})$$

$$\displaystyle\xi_{tk}(\mbox{\boldmath$\Theta$}^{(j-1)})$$

$$\displaystyle\eta_{t-1,l}^{t,k}(\mbox{\boldmath$\Theta$}^{(j-1)})$$

$$\displaystyle\mbox{P}({\bf J}^{\text{\tiny p}}_{t,n}\in V_{k_{n}}|{\bf J}^{\text{\tiny p}}_{t,n^{\prime}}\in V_{k_{n^{\prime}}},\mathcal{Y}_{T},\mbox{\boldmath$\Theta$},n^{\prime}\neq n)$$

$$\displaystyle\zeta_{t,k_{n}}(\mbox{\boldmath$\Theta$}_{\text{s}}^{(j-1)})$$

$$\displaystyle\bigg{|}\mbox{P}({\bf J}_{t,n}^{\text{\tiny p}}\in V_{k_{n}}|{\bf J}_{t,n^{\prime}}^{\text{\tiny p}}\in V_{k_{n}^{\prime}},\mathcal{Y}_{T},\hat{\mbox{\boldmath$\Theta$}}_{\text{s}},n^{\prime}\neq n)-$$

$$\displaystyle\quad\quad\sum_{k_{n^{\prime}}=1}^{K}\mbox{P}({\bf J}^{\text{\tiny p}}_{t,1}\in V_{k_{1}},\ldots,{\bf J}^{\text{\tiny p}}_{t,N}\in V_{k_{N}}|\mathcal{Y}_{T},\hat{\mbox{\boldmath$\Theta$}}_{\text{ns}})\bigg{|}\leq c\varepsilon,$$

Full text

MEG inverse problem,

discrete posterior distribution,

Expectation-Maximization,

spatio-temporal model,

source localization,

arXiv:0000.0000 \startlocaldefs

\endlocaldefs

T1Research is supported by the MOE grants Tier 1 R-155-000-196-114 and Tier 2 R-155-000-184-112 at the National University of Singapore.

, ,

and

1 Introduction

The human brain produces a wide range of bioelectromagnetic signals from various electrical impulses when activated. The signals produced by the neurons in the brain varies from $10$s of femto-Tesla (fT) to $100$s of fT, which is approximately a billion times smaller than the Earth’s magnetic field. Magnetoencephalography (MEG) is a non-invasive imaging technique that is able to detect the weak magnetic fields generated by the neuronal activity within the brain. The MEG recording is able to measure the magnetic fields caused by the neuronal activity inside the brain based on the instrument that is placed close to the scalp. The Superconducting Quantum Inference Devices (SQUIDs) are operated in a magnetically shielded room, and the sensors of SQUIDs are fixed in a one-size-fits-all helmet. During the MEG scanning of SQUIDs, the patient sits under the machine and is restricted from moving. The latest Optically Pumped Magnetometer (OMP) system [4], which is equipped with a customized helmet, allows free and natural movement, including head nodding, stretching, drinking and playing a ball game. With the excellent temporal resolution on a millisecond scale, the MEG has been applied to provide new insights into the neural basis of developmental disorders.

1.1 MEG Inverse Problem

In neuromagnetism, the neuronal current ${\bf J}(r)$ at location $r$ is divided into the primary current ${\bf J}^{\text{\tiny p}}(r)$ and the volume current ${\bf J}^{\text{\tiny v}}(r)$ [7]. Since the primary current ${\bf J}^{\text{\tiny p}}(r)$ is related to the movement of ions and the volume current ${\bf J}^{\text{\tiny v}}(r)$ does not build-up any electric charge, the source of brain activity can be captured by finding the primary current ${\bf J}^{\text{\tiny p}}(r)$. The primary current ${\bf J}^{\text{\tiny p}}(r)$ can then be regarded as current dipoles, and approximated by the summation of $N$ current dipoles,

[TABLE]

where $\delta(\cdot)$ is the Dirac delta function and $Q_{n}$ is a charge dipole at location $r_{n}$, for $n=1,\ldots,N$. The forward problem in neuromagnetism focuses on calculating the magnetic field ${\bf B}$ at location $r$ from a given primary current ${\bf J}^{\text{\tiny p}}(r^{\prime})$ within the brain. Using the quasi-static approximation of the Maxwell’s equations [18], the magnetic field ${\bf B}$ at location $r$ generated by a current dipole $Q_{n}$ is approximated by the Biot-Savart equation,

[TABLE]

where $\mu_{0}$ is the permittivity of free space and $\Omega$ is the brain volume. Thus, the magnetic field at location $r$ of $N$ current dipoles is the summation of ${\bf B}_{n}(r)$ over each dipole $n=1,\ldots,N$.

The MEG inverse problem is to infer the source current given the measured magnetic fields collected from the MEG experiment. However, the general inverse problem is ill-posed. The solution of the inverse problem is not unique since the measured magnetic field could result from an infinite number of possible source currents. This fact makes the MEG inverse problem challenging to solve, and we aim to get a meaningful structure of the source current for the inverse problem. Two types of models have been developed for the MEG inverse problem [3]: equivalent current dipole (ECD) models and distributed source models. The ECD models are based on the assumption that the locations of the current dipoles are unknown and have to be estimated. On the other hand, the distributed source models assume that the measured magnetic fields are generated from the current dipoles with known locations.

1.2 Existing Source Localization Methods

In the literature, two categories of methods focusing on addressing the challenging source localization were proposed. The first category assumes that the current source is static during the MEG scans, which allows us to solve the inverse problem at each time point independently using the quasi-static approximation. Several existing methods were proposed to investigate the current source under the distributed source model and interpret the pattern from the observed magnetic fields. The $L_{2}$ norm and its variations were implemented to solve the distributed source current by using the regularization methods, including the minimum norm estimate [8], minimum current estimate [20], depth-weighted minimum norm estimate [10], and low-resolution electromagnetic tomography algorithm [16]. The MUltiple SIgnal Classification (MUSIC) algorithm [14] is a subspace scanning method, in which the solution is found by searching a single source current through the three-dimensional head volume entirely and projected to an estimate of the signal subspace. The beamforming methods assume that the source currents are uncorrelated, and the goal of the beamformers is to find a set of filter coefficients of the measured magnetic fields, subjected to some constraints. The Linearly Constrained Minimum Variance (LCMV) beamforming method in [22], which was first applied to the inverse problem, is operated by searching a selected region of the head volume to analyze the source current distribution subjected to the minimum variant constraint.

The second category of the methods on source localization incorporates the source current with a time-varying structure. By assuming the variability of source activity, it is able to investigate the source current at each time point $t$, ${\bf J}_{t}^{\text{\tiny p}}$, and provide the temporal evolution of the source current during the MEG scans. The spatio-temporal regularization was utilized in [15] and improves the reconstruction accuracy of the distributed source current. In [13], a dynamic state-space model was proposed to model the movement of the sources, and the Kalman filter (KF) and fixed interval smoother (FIS) were used to solve the high dimensional state estimation. The beamforming method with spatial and temporal effects was proposed to summarize the information of the sources during the voxel-based searching of the head volume, see [25, 26]. Recent work has addressed the source localization of time-varying currents as part of a Bayesian framework. In [2], a Bayesian approach with a maximum a posteriori (MAP) estimator of source activities was built in the distributed source model. The variational Bayesian learning algorithm was derived to reconstruct the distributed sources in the probabilistic generative model, see [19, 6]. In [9], the authors developed a dynamic Maximum a Posterior Expectation-Maximization (dMAP-EM) source localization algorithm based on the KF, FIS, and EM algorithm, to obtain a spatio-temporal distributed solution for the source current. Two sequential important sampling (SIS) [11] based methods, the regular SIS method with rejection and improved SIS method with resampling, were developed in [23] to address the source localization in the ECD models. These authors investigated the source localization by finding the marginal posterior distribution of the source current given the measured magnetic fields, thus providing the variation of the location, orientation, and strength of the source current at each time point.

In neuromagnetism [7], the current dipole $Q_{n}$ is mathematically parameterized with location parameter ${\bf p}_{n}$ and moment parameter ${\bf q}_{n}$. The magnetic fields (1.1), generated from the current dipole $Q_{n}$, can be approximated by the Biot-Savart in a discrete matrix form,

[TABLE]

where ${\bf L}_{n}(r,{\bf p}_{n})$ is the lead field. In the distributed source models, the location parameter ${\bf p}_{n}$ of the source current is assumed to be known, thus the lead field in (1.2) can be calculated from the forward model. In this case, the magnetic fields linearly depend on the estimated moment parameter ${\bf q}_{n}$. Several approaches, such as the FK and FIS, were proposed to estimate the moment parameter ${\bf q}_{n}$, and they work well under the linear dynamic system. However, the magnetic fields also non-linearly depend on the unknown location parameter ${\bf p}_{n}$, thus the recovery of the location parameter usually involves non-linear optimization, in which applying the Kalman filter would degrade the performance [1]. The existing work on time varying source current with unknown location parameter is limited, and this motivates us to develop new approaches to investigate the source distribution for the location parameter.

1.3 Goal of this Paper

The goal of this paper is to invent a new Bayesian framework to find the posterior distribution of the source current ${\bf J}^{\text{\tiny p}}_{t}$ at time point $t$, given the entire collection of measurements $\mathcal{Y}_{T}$, which consists of measurements ${\bf Y}_{t}$ for $1\leq t\leq T$. The posterior distribution $p({\bf J}^{\text{\tiny p}}_{t}|\mathcal{Y}_{T})$ for the source current ${\bf J}^{\text{\tiny p}}_{t}$ can be interpreted as a solution for the MEG inverse problem. In contrast to the existing literature, our proposed methodology is novel based on the following two aspects. First, we develop a discrete approach for computing the discrete posterior distribution of the source current ${\bf J}^{\text{\tiny p}}_{t}$, and the discrete posterior distribution is used to approximate the continuous posterior distribution $p({\bf J}^{\text{\tiny p}}_{t}|\mathcal{Y}_{T})$. The SIS schemes in [23] investigated the source distribution by numerically sampling the continuous posterior distribution. However, the posterior distribution in [23] does not have an analytically tractable closed form, and the sampling procedure may suffer from slow convergence due to the high dimensionality of the measurements. In comparison with the sampling schemes, our method gives the discrete posterior distribution for the source current with a closed form, and are able to calculate it directly even when the dimensionality of the measurements is large. Second, we use the entire collection of measurements to investigate the source distribution $p({\bf J}^{\text{\tiny p}}_{t}|\mathcal{Y}_{T})$ instead of the source distribution $p({\bf J}^{\text{\tiny p}}_{t}|\mathcal{Y}_{t})$, using only the measurements up to the time point that we are interested in. The MEG allows for a real-time recording of the brain activity on a millisecond scale. For each time point $t$, the past measurements and the future measurements both reflect the trajectory of the time-varying source ${\bf J}^{\text{\tiny p}}_{t}$, for $1\leq t\leq T$. In contrast to previous related approaches [2, 19, 23], we utilize the entire set of measurements to recover the location of the source current.

For the proposed discrete approach, we focus on the selected three-dimensional region of interest (ROI), and the ROI is subsequently discretized into $K$ voxels $\{V_{k}\}_{k=1}^{K}$. Then, we calculate the discrete posterior distribution $\mbox{P}({\bf J}_{t}^{\text{\tiny p}}\in V_{k}|\mathcal{Y}_{T})$ of the source current ${\bf J}_{t}^{\text{\tiny p}}$ at each time point $t$, for all $1\leq k\leq K$. Figure 7 presents the posterior distribution for the location parameter of a single source at six selected time points, where the probabilities are highlighted in different colors. The discrete posterior distribution indicates that the source would appear in the voxel with a corresponding probability. The region with non-zero probabilities can be interpreted as the activated area of the source current. Thus, the discrete approach is able to provide the source distribution on time evolution during the MEG scanning.

In order to calculate the discrete posterior distribution, the EM algorithm with incomplete data is implemented to estimate the parameters in the source model. We further develop a switch procedure and dynamic procedure to implement the proposed discrete approach. The switch procedure is proposed to deal with the case involving multiple sources, and the dynamic procedure is developed to balance the estimation accuracy and computational efficiency when no available information on the ROI is provided. With the proposed approach, we will be able to investigate both MEG and EEG recordings that contain valuable time-sensitive information and shed light on the time evolution of the source localization.

1.4 Outline of this Paper

In Section 2, we present the methodology of the discrete approach used to recover the source distribution. First, we utilize a dynamic spatio-temporal model to reformulate the source localization problem in Section 2.1. After which, we propose the discrete approach for the calculation of the posterior distribution for the location parameter of the source current in Section 2.2. An estimation procedure for the parameters, which are introduced by the source model, will be presented in Section 2.3. We further develop the switch procedure and dynamic procedure to implement the discrete approach for the calculation of the posterior distribution in Section 2.4 and Section 2.5. Simulation studies are described in Section 3. In Section 4, a MEG data application is presented. An extension to the EEG data is illustrated in Section 5. A short discussion and concluding remarks are given in Section 6.

2 Methodology

2.1 A Dynamic Spatio-Temporal Model

In an MEG experiment, the observed magnetic fields are scanned at $L$ sensors, and the data is recorded for a fixed time period, $T$ milliseconds. Let ${\bf Y}_{t}=(Y_{t,1},\ldots,Y_{t,L})^{\rm T}$ be the measurements of $L$ sensors collected at time point $t$, $1\leq t\leq T$, and $\mathcal{Y}_{T}=({\bf Y}_{1},\ldots,{\bf Y}_{T})$ be the entire collection of measurements of the experiment. The magnetic field measured from the $l$-th sensor at time point $t$ is

[TABLE]

where $U_{t,l}$ denotes the noise of the measurements. In this paper, we consider a horizontally layered conductor, and the magnetic field is sensitive only to the tangential component of the source. Thus, the magnetic field ${\bf B}_{l}({\bf J}^{\text{\tiny p}}_{t})$ in (1.1) is defined as

[TABLE]

where ${\bf J}^{\text{\tiny p}}_{t}=({\bf p}_{t}^{\rm T},{\bf q}_{t}^{\rm T})^{\rm T}$, ${\bf p}_{t}=(p_{t,1},p_{t,2},p_{t,3})^{\rm T}$ contains the location parameters, ${\bf q}_{t}=(q_{t,1},q_{t,2},q_{t,3})^{\rm T}$ contains the moments and strength of the source at time $t$, ${\bf r}_{l}$ is the location of the $l$-th sensor, and ${\bf e}_{z}=(0,0,1)^{\rm T}$.

Remark 2.1**.**

All the methodology developed in this paper is legitimately extendable to the EEG source analysis. Speaking of the extension to EEG, the potential field generated by the source currents is used to replace the magnetic field defined in (2.1). Similarly, the potential field generated by a source ${\bf J}^{\text{\tiny p}}_{t}$ at the $l$-th sensor is given by,

[TABLE]

where $\sigma$ is the conductivity and $1\leq l\leq L$.

Let ${\bf B}({\bf J}^{\text{\tiny p}}_{t})=(B_{1}({\bf J}^{\text{\tiny p}}_{t}),\ldots,B_{L}({\bf J}^{\text{\tiny p}}_{t}))^{\rm T}$ be the collection of magnetic fields with source current ${\bf J}^{\text{\tiny p}}_{t}$ generated at $L$ sensors. Then, the general framework of the MEG measurement ${\bf Y}_{t}$ is given by

[TABLE]

where ${\bf U}_{t}=(U_{t,1},\dots,U_{t,L})^{\rm T}$ contains the noise of the measurements. We assume that the noises $U_{t,l}$, $1\leq l\leq L$, are uncorrelated between every pair of sensors and homogeneous. For simplicity, we assume that the noise is Gaussian, that is, ${\bf U}_{t}\sim\mathcal{N}({\bf 0},{\bf V})$, where ${\bf V}=\text{diag}(\sigma^{2},\ldots,\sigma^{2})$. To investigate the source distribution, we consider a dynamic spatio-temporal model for the source ${\bf J}^{\text{\tiny p}}_{t}$. The first order auto-regression model is given by

[TABLE]

where {\bf J}^{\text{\tiny p}}_{1}\sim\mathcal{N}(\mbox{\boldmath\mu}_{0},\mbox{\boldmath\Sigma}_{0}), and the evolution noise ${\bf Z}_{t}$ is assumed to be Gaussian, {\bf Z}_{t}\sim\mathcal{N}({\bf 0},\mbox{\boldmath\Sigma}), \mbox{\boldmath\Sigma}=\text{diag}(\sigma_{1}^{2},\ldots,\sigma_{6}^{2}). Let $\mathcal{J}^{\text{\tiny p}}_{T}$ be the collection of source ${\bf J}^{\text{\tiny p}}_{t}$ for $1\leq t\leq T$. The causal relationship between the MEG measurements $\mathcal{Y}_{T}$ and source $\mathcal{J}^{\text{\tiny p}}_{T}$ is described in Figure 1.

Throughout the paper, the framework is mainly based on a single source and can be generalized to the case with multiple sources. We extend ${\bf J}^{\text{\tiny p}}_{t}$ to be the collection of $N$ current sources, where ${\bf J}^{\text{\tiny p}}_{t}=({{\bf J}^{\text{\tiny p}}_{t,1}}^{\rm T},\ldots,{{\bf J}^{\text{\tiny p}}_{t,N}}^{\rm T})^{\rm T}$, and ${\bf J}^{\text{\tiny p}}_{t,n}$ is the $n$-th source parameterized with $({\bf p}_{t,n}^{\rm T},{\bf q}_{t,n}^{\rm T})^{\rm T}$ for $1\leq n\leq N$. Then, the magnetic field in (2.1) is generated from $N$ sources and is given by ${\bf B}_{l}({\bf J}^{\text{\tiny p}}_{t})=\sum_{n=1}^{N}{\bf B}_{l}({\bf J}^{\text{\tiny p}}_{t,n})$. We further assume that the $N$ sources are uncorrelated, and each source ${\bf J}^{\text{\tiny p}}_{t,n}$ is modeled with the first order auto-regression model,

[TABLE]

where {\bf J}^{\text{\tiny p}}_{1,n}\sim\mathcal{N}(\mbox{\boldmath\mu}_{0,n},\mbox{\boldmath\Sigma}_{0,n}), {\bf Z}_{t,n}\sim\mathcal{N}({\bf 0},\mbox{\boldmath\Sigma}_{n}). Thus, the general framework of $N$ uncorrelated sources ${\bf J}^{\text{\tiny p}}_{t}$ is consistent with model (2.4), where the parameters \mbox{\boldmath\mu}_{0}, \mbox{\boldmath\Sigma}_{0}, $\Sigma$, $A$ and ${\bf b}$ contain the information of $N$ sources correspondingly. To be precise, $A=\text{diag}(A_{1},\ldots,A_{N})$, ${\bf b}=({\bf b}_{1}^{\rm T},\ldots,{\bf b}_{N}^{\rm T})^{\rm T}$, \mbox{\boldmath\mu}_{0}=(\mbox{\boldmath\mu}_{0,1}^{\rm T},\ldots,\mbox{\boldmath\mu}_{0,N}^{\rm T})^{\rm T}, \mbox{\boldmath\Sigma}_{0}=\text{diag}(\mbox{\boldmath\Sigma}_{0,1},\ldots,\mbox{\boldmath\Sigma}_{0,N}), and \mbox{\boldmath\Sigma}=\text{diag}(\mbox{\boldmath\Sigma}_{1},\ldots,\mbox{\boldmath\Sigma}_{N}).

Let \mbox{\boldmath\Theta}=\{\mbox{\boldmath\mu}_{0},\mbox{\boldmath\Sigma}_{0},A,{\bf b},\mbox{\boldmath\Sigma},{\bf V}\} be the list of parameters introduced in the framework of measurements (2.3) and source model (2.4). From the causal relationship described in Figure 1, we note that the sequences of source $\mathcal{J}^{\text{\tiny p}}_{T}$ and measurements $\mathcal{Y}_{T}$ have the following Markov properties.

p({\bf J}^{\text{\tiny p}}_{t}|\mathcal{J}^{\text{\tiny p}}_{t-1},\mbox{\boldmath\Theta})=p({\bf J}^{\text{\tiny p}}_{t}|{\bf J}^{\text{\tiny p}}_{t-1},\mbox{\boldmath\Theta}), for $2\leq t\leq T$. 2. 2.

p({\bf J}^{\text{\tiny p}}_{t}|{\bf J}^{\text{\tiny p}}_{t-1},\mathcal{Y}_{t-1},\mbox{\boldmath\Theta})=p({\bf J}^{\text{\tiny p}}_{t}|{\bf J}^{\text{\tiny p}}_{t-1},\mbox{\boldmath\Theta}), for $2\leq t\leq T$. 3. 3.

p(\mathcal{Y}_{t}|\mathcal{J}^{\text{\tiny p}}_{T},\mbox{\boldmath\Theta})=p({\bf Y}_{t}|{\bf J}^{\text{\tiny p}}_{t},\mbox{\boldmath\Theta}), for $1\leq t\leq T$.

Throughout the paper, we use $p$ as a genetic symbol for continuous probability distribution.

2.2 A Discrete Approach

Given the entire magnetic measurement $\mathcal{Y}_{T}$, we aim to investigate the source distribution for location parameter ${\bf p}_{t}$ by computing the posterior distribution of ${\bf J}^{\text{\tiny p}}_{t}$, for $1\leq t\leq T$. From Bayes theorem, we have

[TABLE]

where $\mathcal{Y}_{t}$ is the collection of the measurements up to time point $t$ and $\mathcal{Y}_{T\setminus t}$ contains the remaining measurements from time point $t+1$ to $T$.

In this section, we build a discrete approach to approximate the posterior distribution in (2.5), reformulate the discrete posterior distribution with a closed form, and present the forward-backward algorithm to further compute it. The location parameter of interest is ${\bf p}_{t}$ of the source current ${\bf J}_{t}^{\text{\tiny p}}$, for $1\leq t\leq T$, thus the moment and strength parameter ${\bf q}_{t}$ is fixed for all the time points throughout this paper.

With an available ROI, the movement of location parameter ${\bf p}_{t}$ of source ${\bf J}^{\text{\tiny p}}_{t}$ is assumed to be restricted within the ROI at all times. For the discrete approach, each dimension of the ROI is discretized with mesh grid $K_{i}$, $i=1,2,3$. The mesh grids are used to construct a sequence of voxels $\{V_{k}\}_{k=1}^{K}$ and approximate the three-dimensional ROI, where $V_{k}$ is the $k$-th voxel with its center ${\bf c}_{k}$ and $K=K_{1}\cdot K_{2}\cdot K_{3}$. To investigate the source current ${\bf J}^{\text{\tiny p}}_{t}$, we introduce a corresponding set of binary indicator variables $v_{tk}\in\{0,1\}$, where $k=1,\ldots,K$, describing which of the $K$ voxels the source ${\bf J}^{\text{\tiny p}}_{t}$ is located in, so that if ${\bf J}^{\text{\tiny p}}_{t}\in V_{k}$ then $v_{tk}=1$, and $v_{tk^{\prime}}=0$ for $k^{\prime}\neq k$. Applying the coding scheme, we can then approximate the continuous probabilities introduced in (2.3) and (2.4) using the corresponding discrete probabilities. For the continuous probability p({\bf J}^{\text{\tiny p}}_{1}|\mbox{\boldmath\Theta}), there exists $1\leq k^{\ast}\leq K$, such that

[TABLE]

Similarly, we have the following approximations

[TABLE]

for $2\leq t\leq T$, and

[TABLE]

for $1\leq t\leq T$.

Under the discrete approach, the posterior distribution in (2.5) can be approximated by

[TABLE]

for $1\leq k\leq K$ and $1\leq t\leq T$. Let \alpha_{tk}(\mbox{\boldmath\Theta})=\mbox{P}(v_{tk}=1,\mathcal{Y}_{t}|\mbox{\boldmath\Theta})/p(\mathcal{Y}_{t}|\mbox{\boldmath\Theta}), \beta_{tk}(\mbox{\boldmath\Theta})=\mbox{P}(\mathcal{Y}_{T\setminus t}|v_{tk}=1,\mbox{\boldmath\Theta})/p(\mathcal{Y}_{T\setminus t}|\mbox{\boldmath\Theta}), and we have

[TABLE]

Therefore, the calculation of the discrete posterior distribution (2.9) consists of two parts, a filtering procedure on \alpha_{tk}(\mbox{\boldmath\Theta}) using the past measurements $\mathcal{Y}_{t}$ and a smoothing procedure on \alpha_{tk}(\mbox{\boldmath\Theta}) using the remaining measurements $\mathcal{Y}_{T\setminus t}$.

The forward-backward algorithm [17] is an efficient inference algorithm, which computes the posterior distribution of the hidden state variables given the entire set of measurements in two passes. The forward-backward algorithm for computing the discrete posterior distribution in (2.9) is summarized in Table 1 and Table 2. In Table 1, we start the forward recursion from the first time point $t=1$ and compute the filtering posterior distribution \alpha_{1k}(\mbox{\boldmath\Theta}) given $\Theta$, for $1\leq k\leq K$. After which, we compute the filtering posterior distribution \alpha_{tk}(\mbox{\boldmath\Theta}) using the previous filtering posterior distribution \alpha_{t-1,l}(\mbox{\boldmath\Theta}), for $t=2,\ldots,T$. Table 2 illustrates the backward procedure to calculate the smoothing posterior distribution \beta_{tk}(\mbox{\boldmath\Theta}). We start the calculation from the last time point $t=T$ and initialize \beta_{Tk}(\mbox{\boldmath\Theta})=1 for $1\leq k\leq K$. Then, we calculate the smoothing posterior distribution \beta_{tk}(\mbox{\boldmath\Theta}) using the smoothing posterior distribution \beta_{t+1,l}(\mbox{\boldmath\Theta}), for $t=T-1,\ldots,1$. From the output of the forward-backward algorithm, the discrete posterior distribution of the source current ${\bf J}^{\text{\tiny p}}_{t}$ is given in (2.9).

2.3 Parameter Estimation

Under the models (2.3) and (2.4), the posterior distribution defined in (2.5) depends on the parameter $\Theta$. In this section, we apply the EM algorithm [5] to find a MLE \hat{\mbox{\boldmath\Theta}} of the parameter $\Theta$ with incomplete MEG data $(\mathcal{Y}_{T},\mathcal{J}^{\text{\tiny p}}_{T})$, as we have no access to the collection of the source $\mathcal{J}^{\text{\tiny p}}_{T}$ during the MEG scans. The optimization problem is defined as

[TABLE]

where \ell(\mbox{\boldmath\Theta},\mathcal{Y}_{T}) is the log likelihood function of parameter $\Theta$, given the entire set of measurements $\mathcal{Y}_{T}$. Under the discrete approach in Section 2.2, the unobserved source $\mathcal{J}^{\text{\tiny p}}_{T}$ is assumed to be the discrete variable, thus the log likelihood function in (2.10) can be written into the following form

[TABLE]

where

[TABLE]

under the Markov properties of $\mathcal{Y}_{T}$ and $\mathcal{J}^{\text{\tiny p}}_{T}$. Applying Jensen’s inequality, we have

[TABLE]

where $q(\mathcal{J}^{\text{\tiny p}}_{T})$ is a probability distribution on the unobserved variables $\mathcal{J}^{\text{\tiny p}}_{T}$. From (2.3), we have

[TABLE]

for any probability distribution $q(\mathcal{J}^{\text{\tiny p}}_{T})$. When q(\mathcal{J}^{\text{\tiny p}}_{T})=\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath\Theta}), the equality in (2.3) holds,

[TABLE]

To maximize the log-likelihood function \ell(\mbox{\boldmath\Theta},\mathcal{Y}_{T}) in (2.10), we construct an alternating EM algorithm, maximizing L(\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath\Theta}),\mbox{\boldmath\Theta},\mathcal{Y}_{T}) defined in (2.13).

We start the EM algorithm from a reasonable initialization \mbox{\boldmath\Theta}^{(0)}. Let \mathcal{L}(\mbox{\boldmath\Theta}|\mbox{\boldmath\Theta}^{(j-1)}):=L(\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath\Theta}^{(j-1)}),\mbox{\boldmath\Theta},\mathcal{Y}_{T}), $j=1,2,\ldots$. For the following iterations, we compute the posterior probability \mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath\Theta}^{(J-1}) to maximize \mathcal{L}(\mbox{\boldmath\Theta}|\mbox{\boldmath\Theta}^{(j-1)}) in the E-step. It is noted that we cannot treat the summation over $\mathcal{J}_{T}^{\text{\tiny p}}$ in (2.13) at each time point independently. Since the expectation of the binary variable $v_{tk}$ is just the probability that it takes the value $1$, we have

[TABLE]

for $1\leq t\leq T$, $1\leq k\leq K$, and

[TABLE]

for $2\leq t\leq T$, $1\leq k\leq K$. For the first term of \mathcal{L}(\mbox{\boldmath\Theta}|\mbox{\boldmath\Theta}^{(j-1)}) defined in (2.13), we have

[TABLE]

where the complete likelihood function p(\mathcal{Y}_{T},\mathcal{J}^{\text{\tiny p}}_{T}|\mbox{\boldmath\Theta}) is approximated by the discrete distributions (2.6)-(2.8). For the $j$-th iteration, we let \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}):=\mbox{P}(v_{tk}=1|\mathcal{Y}_{T},\mbox{\boldmath\Theta}^{(j-1)}) be the intermediate discrete posterior distribution, and \eta_{t-1,l}^{t,k}(\mbox{\boldmath\Theta}^{(j-1)}):=\mbox{P}(v_{t-1,l}=1,v_{tk}=1| \mathcal{Y}_{T},\mbox{\boldmath\Theta}^{(j-1)}) be the intermediate discrete joint posterior distribution. Instead of targeting the posterior distribution $\mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T}$, \mbox{\boldmath\Theta}^{(j-1)}), we compute the intermediate posterior distributions \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) and \eta_{t-1,l}^{t,k}(\mbox{\boldmath\Theta}^{(j-1)}) during the E-step at the $j$-th iteration. To be precise, we obtain \alpha_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) and \beta_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) from the forward-backward algorithm given the intermediate estimate \mbox{\boldmath\Theta}^{(j-1)}. Then, we have

[TABLE]

for $1\leq t\leq T$, and

[TABLE]

where $2\leq t\leq T$, $1\leq k\leq K$.

In the M-step of the $j$-th iteration, we update \mbox{\boldmath\Theta}^{(j)}, in which the intermediate posterior distributions \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) and \eta_{t-1,l}^{t,k}(\mbox{\boldmath\Theta}^{(j-1)}) are treated as constant. The maximization of \mathcal{L}(\mbox{\boldmath\Theta}|\mbox{\boldmath\Theta}^{(j-1)}) is equivalent to maximize (2.16), and it gives the closed form of the updates \mbox{\boldmath\Theta}^{(j)}, see Table 3.

We perform the E-step and M-step until the function \mathcal{L}(\mbox{\boldmath\Theta}|\mbox{\boldmath\Theta}^{(j-1)}) defined in (2.13) converges. The EM algorithm is summarized in Table 4. Under the following regularity conditions,

1. (C1).

The parameter space $\Omega$ is an open set in the Euclidean space. 2. (C2).

The variables ${\bf Y}_{t}$, $1\leq t\leq T$ are independent and identically distributed with density p(y|\mbox{\boldmath\Theta}). There exists positive constants $c$ and $C$, such that c\leq p(y|\mbox{\boldmath\Theta})\leq C. Moreover, the density function p(y|\mbox{\boldmath\Theta}) is differentiable with continuous first derivative. 3. (C3).

The level set \Omega_{\mbox{\boldmath\Theta}}=\{\mbox{\boldmath\Theta}^{\prime}\in\Omega:\ell(\mbox{\boldmath\Theta}^{\prime},\mathcal{Y}_{T})\geq\ell(\mbox{\boldmath\Theta},\mathcal{Y}_{T})\} is compact. 4. (C4).

The distribution of the unobserved variables \mbox{P}(\mathcal{J}^{\text{\tiny p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath\Theta}) has the same support for all \mbox{\boldmath\Theta}\in\Omega. 5. (C5).

The function Q(\mbox{\boldmath\Theta}^{\prime}|\mbox{\boldmath\Theta}) is continuous in both $\Theta$ and \mbox{\boldmath\Theta}^{\prime} and differentiable in \mbox{\boldmath\Theta}^{\prime}. 6. (C6).

All the stationary points in $\mathcal{S}_{\ell}$ are isolated, where $\mathcal{S}_{\ell}$ denotes the stationary points of the log likelihood function \ell(\mbox{\boldmath\Theta},\mathcal{Y}_{T}). 7. (C7).

For all \mbox{\boldmath\Theta}\in\mathcal{S}_{\ell} there exists a unique global maximum of Q(\cdot|\mbox{\boldmath\Theta}),

we prove the convergence of the EM sequence \{\mbox{\boldmath\Theta}^{(j)}\} in Theorem 2.1 and the proof can be found in Appendix A.

Theorem 2.1**.**

The iterative procedure of EM algorithm does not cause a decrease in the log likelihood function \ell(\mbox{\boldmath\Theta},\mathcal{Y}_{T}). Furthermore, we assume that (C1) - (C7) hold. Then, for any starting value \{\mbox{\boldmath\Theta}^{(0)}\}, the EM sequence \{\mbox{\boldmath\Theta}^{(j)}\}, \mbox{\boldmath\Theta}^{(j)}\to\mbox{\boldmath\Theta}^{\ast}, when $j\to\infty$, for some stationary point \mbox{\boldmath\Theta}^{\ast}\in\mathcal{S}_{\ell}.

2.4 A Switch Procedure

When we apply the discrete approach to the case with multiple sources, the possible states for $N$ sources go to $K^{N}$, where $K$ is the number of discrete voxels. In the Supplementary Materials, we compare the possible states for $N$ sources with mesh grids in the discretization. Even with the mesh grids $K_{i}=6$, $i=1,2,3$, the calculation for the discrete posterior distribution needs to include more than $40000$ possible states, which makes the computation procedure impossible. Thus, a method that enables a more achievable posterior calculation of multiple sources is desirable.

In this section, we propose a switch procedure to reduce the computational complexity for the case with multiple sources. Applying the Bayes rule, we have

[TABLE]

for $1\leq n\leq N$. To reduce the computational complexity, the marginal posterior probability \mbox{P}({\bf J}^{\text{\tiny p}}_{t,n}\in V_{k_{n}}|{\bf J}^{\text{\tiny p}}_{t,n^{\prime}}\in V_{k_{n^{\prime}}},\mathcal{Y}_{T},\mbox{\boldmath\Theta},n^{\prime}\neq n), $1\leq n\leq N$, is used to approximate the posterior distribution of $N$ sources, and we calculate the marginal posterior distribution of the source ${\bf J}^{\text{\tiny p}}_{t,n}$ by assuming that the states of other sources ${\bf J}^{\text{\tiny p}}_{t,n^{\prime}}$ are known. During the EM iterations, we let \zeta_{t,k_{n}}(\mbox{\boldmath\Theta}_{\text{s}}^{(j-1)})=\mbox{P}(v_{t,k_{n}}=1|{\bf J}^{\text{\tiny p}}_{t,n^{\prime}}\in V_{k_{n^{\prime}}},\mathcal{Y}_{T},\mbox{\boldmath\Theta}_{\text{s}}^{(j-1)},n^{\prime}\neq n) at the $j$-th iteration, where \{\mbox{\boldmath\Theta}_{\text{s}}^{(j)}\} is the EM sequence with the switch procedure. For $1\leq n\leq N$, we compute the marginal posterior probability by

[TABLE]

When the EM iteration converges, we will use the estimated parameter \hat{\mbox{\boldmath\Theta}}_{\text{s}} to calculate the posterior distribution \mbox{P}({\bf J}^{\text{\tiny p}}_{t,1}\in V_{k_{1}},\ldots,{\bf J}^{\text{\tiny p}}_{t,N}\in V_{k_{N}}|\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{s}}), which is approximated by the marginal posterior distribution \mbox{P}({\bf J}^{\text{\tiny p}}_{t,n}\in V_{k_{n}}|{\bf J}^{\text{\tiny p}}_{t,n^{\prime}}\in V_{k_{n^{\prime}}},\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{s}},n^{\prime}\neq n) in the switch procedure.

The posterior distribution from the switch procedure is compared with the one from the non-switch procedure in Theorem 2.2 and the proof can be found in Appendix B.

Theorem 2.2**.**

Under regularity conditions (C1) and (C7), there exists $\varepsilon>0$, such that

[TABLE]

where $1\leq n\leq N$, $c$ is a positive constant, \mbox{P}({\bf J}_{t,n}^{\text{\tiny p}}\in V_{k_{n}}|{\bf J}_{t,n^{\prime}}^{\text{\tiny p}}\in V_{k_{n}^{\prime}},\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{s}},n^{\prime}\neq n) is the marginal posterior distribution of source ${\bf J}_{t,n}^{\text{\tiny p}}$ obtained from the switch procedure, \hat{\mbox{\boldmath\Theta}}_{\text{ns}} is the estimate obtained from the non-switch procedure, and \mbox{P}({\bf J}^{\text{\tiny p}}_{t,1}\in V_{k_{1}},\ldots,{\bf J}^{\text{\tiny p}}_{t,N}\in V_{k_{N}}|\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{ns}}) is the posterior distribution of $N$ sources from the non-switch procedure.

2.5 A Dynamic Procedure

When no information on the ROI is available, we have to discretize the whole head model. Implementing the discrete approach to the parameter estimation procedure, we note that the estimate depends on the discretization. Increasing mesh girds will improve the estimation accuracy of calculating the discrete posterior distribution of the source during the EM iterations. However, it will also increase the computational complexity. This motivates us to develop a dynamic procedure to implement the discrete approach in order to balance the estimation accuracy and computational complexity when no information on the ROI is available.

We assume that a true ROI, $\text{ROI}_{0}$, exists and is assumed to restricted within the head model. The movement of the source current ${\bf J}_{t}^{\text{\tiny p}}$ is assumed to be restricted within $\text{ROI}_{0}$ at all times. In the EM algorithm with the discrete approach, calculating the discrete posterior distributions \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) and \eta_{t-1,l}^{tk}(\mbox{\boldmath\Theta}^{(j-1)}) depends on the discretization in the E-step of the $j$-th iteration. For the dynamic procedure, we assume that a shrunken ROI, $\text{ROI}^{(j)}$, exists and covers the true ROI. Utilizing the shrunken ROI, $\text{ROI}^{(j)}$, with increased mesh grids, the dynamic procedure will improve the estimation accuracy of calculating the discrete posterior distributions \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) and \eta_{t-1,l}^{tk}(\mbox{\boldmath\Theta}^{(j-1)}) in the E-step of the $j$-th iteration. Then, the posterior distributions are used to update \mbox{\boldmath\Theta}^{(j)} and the update will also be improved in the following M-step. By implementing the dynamic procedure, we expect to obtain an accurate estimate \hat{\mbox{\boldmath\Theta}}_{\text{d}} with a shrunken ROI.

In the dynamic procedure, we calculate the intermediate posterior distribution \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) for the location parameter ${\bf p}_{t}$ based on the current ROI in the $j$-th iteration. After which, we use the intermediate marginal posterior distribution \xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)}) to obtain a shrunken ROI for the $(j+1)$-th iteration. This new ROI is constructed by $\text{ROI}^{(j+1)}=I_{1}^{(j+1)}\times I_{2}^{(j+1)}\times I_{3}^{(j+1)}$, where the one-dimensional interval is given by

[TABLE]

where \mu^{(j)}_{t,i}=\sum_{k=1}^{K}\xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)})c_{k,i}, \sigma^{(j)}_{t,i}=\text{sqrt}(\sum_{k=1}^{K}\xi_{tk}(\mbox{\boldmath\Theta}^{(j-1)})(c_{k,i}-\mu^{(j)}_{t,i})^{2}), $c_{k,i}$ is the $i$-th component of ${\bf c}_{k}$, and $i=1,2,3$. With the shrunken ROI, we also increase the mesh grids $K_{i}$, $i=1,2,3$, during the iterative procedures. When the dynamic EM algorithm converges, an estimate \hat{\mbox{\boldmath\Theta}}_{\text{d}} will be available with a shrunken ROI.

3 Simulation Study

3.1 MEG Data Generation

In the simulation study, we considered a single sphere head model (centered at the origin with radius $10$ cm), and simulated $102$ magnetometers which are randomly placed on the upper part of the head. Movement of the current sources inside the head model was restricted. In order to focus on the location parameters of the sources, we fixed some parameters (initial distribution parameters \mbox{\boldmath\mu}_{0}, \mbox{\boldmath\Sigma}_{0}, and noise parameters $\Sigma$, ${\bf V}$) in the model (2.3) and (2.4). The total length of each simulation was $100$ time points, and of interest at any time point is the discrete posterior distribution of the sources.

3.2 Simulated Case $1$

In this example, we consider a single source which moves in $3$ dimensions $(x,y,z)$ in the head model. The parameters of the simulated source are summarized in Table 5.

In the simulation, we applied the EM algorithm with the dynamic procedure to estimate the parameter \mbox{\boldmath\Theta}=\{A,{\bf b}\} in the single source model. We chose the upper head model as the initial ROI, and let the initial mesh grids $K_{i}^{(1)}=10$, $i=1,2,3$. Since we applied the discrete approach to the simulated data, we started the dynamic EM algorithm by discretizing the initial ROI with the initial mesh grids. In the $j$-th iteration, the mesh grid was increased by $K_{i}^{(j+1)}=K_{i}^{(j)}+1$, $i=1,2,3$, and the ROI was shrunk for the following iteration. The iteration procedure was performed until the conditional expectation in (2.16) converged, and we obtained the MLEs \hat{\mbox{\boldmath\Theta}}_{\text{d}}=\{\hat{A}_{\text{d}},\hat{{\bf b}}_{\text{d}}\}, where

$\hat{A}_{\text{d}}=\begin{pmatrix}\hat{A}_{\text{d}}^{\ast}&{\bf 0}_{3\times 3}\\ {\bf 0}_{3\times 3}&{\bf 0}_{3\times 3}\end{pmatrix}$, $\hat{A}_{\text{d}}^{\ast}=\begin{pmatrix}0.6969&0.0260&-0.0024\\ 0.0479&0.8352&0.0382\\ -0.0207&0.0029&0.9035\end{pmatrix}$

, and $\hat{{\bf b}}_{\text{d}}=(0.9173,-0.6657,0.3194,3.0032,2.9536,3.0488)^{\rm T}$. When we obtained the estimated parameter, we also obtained a shrunken ROI $[-1.7707,3.1606]$ cm $\times[-3.0992,1.2701]$ cm $\times[1.4235,4.7374]$ cm. Utilizing the estimated parameter and the shrunken ROI, we computed the discrete posterior distribution \mbox{P}({\bf J}^{\text{\tiny p}}_{t}|\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{d}}). To visualize the discrete posterior distribution, we plot the marginal posterior means for location parameter ${\bf p}_{t}$ of the source in Figure 2.

We also compared the numerical result from the dynamic procedure with the result from the non-dynamic EM procedure. For the non-dynamic EM algorithm, we implemented the discrete approach by discretizing the initial ROI with initial mesh grids throughout the simulations. In this case, we obtained the estimates

$\hat{A}_{\text{nd}}=\begin{pmatrix}\hat{A}_{\text{nd}}^{\ast}&{\bf 0}_{3\times 3}\\ {\bf 0}_{3\times 3}&{\bf 0}_{3\times 3}\end{pmatrix}$, $\hat{A}^{\ast}_{\text{nd}}=\begin{pmatrix}0.7280&-0.0018&0.0052\\ -0.0019&0.7888&0.0208\\ -0.0689&-0.0258&0.8988\end{pmatrix}$

, and $\hat{{\bf b}}_{\text{nd}}=(0.7980,-0.6292,0.4197,3.0032$, $2.9536,3.0488)^{\rm T}$. The posterior means for the location parameter of the source are plotted in Figure 3. To compare the estimation accuracy from the two procedures, we repeated the simulation four times and computed the estimation error in maximum norm. From Table 6, we note that the estimates from the dynamic EM procedure are more accurate than the non-dynamic EM procedure.

3.3 Simulated Case $2$

In addition to Case 1, a case of two sources was performed. In the simulation, the two sources were assumed to be uncorrelated and were allowed to move in $3$ dimensions $(x,y,z)$ in the head model. The parameters of the simulated sources are summarized in Table 7.

To deal with the simulated data with two sources, the dynamic procedure and switch procedure were applied to implement the EM algorithm for estimating the parameters \mbox{\boldmath\Theta}=\{A,{\bf b}\} in the source model. When the EM algorithm converged, we obtained the estimates $\hat{A}=\begin{pmatrix}\hat{A}_{1}&{\bf 0}_{6\times 6}\\ {\bf 0}_{6\times 6}&\hat{A}_{2}\end{pmatrix}$, and $\hat{{\bf b}}=(\hat{{\bf b}}_{1}^{\rm T},\hat{{\bf b}}_{2}^{\rm T})^{\rm T}$, where $\hat{A}_{i}=\begin{pmatrix}\hat{A}_{i}^{\ast}&{\bf 0}_{3\times 3}\\ {\bf 0}_{3\times 3}&{\bf 0}_{3\times 3}\end{pmatrix}$, $i=1,2$,

$\hat{A}_{1}^{\ast}=\begin{pmatrix}0.2660&0.0942&-0.1494\\ 0.0283&0.5690&0.2668\\ 0.0884&0.0570&0.8590\end{pmatrix}$

,

$\hat{A}_{2}^{\ast}=\begin{pmatrix}0.2501&-0.0450&0.0243\\ -0.0953&0.7234&0.0916\\ -0.0481&0.0255&0.8550\end{pmatrix}$

,

$\hat{{\bf b}}_{1}=(3.6058,-2.5997,0.2934,2.9535,3.0124,2.9188)^{\rm T}$

, and

$\hat{{\bf b}}_{2}=(2.5589,-0.9985,0.6659,3.0527$, $3.0327,3.0541)^{\rm T}$

. The discrete posterior distributions of the sources were calculated using the estimated parameter \hat{\mbox{\boldmath\Theta}}, and the posterior means for the location parameter of the two sources are plotted in Figure 4.

To compare the result with the non-dynamic procedure, we performed the simulation with the non-dynamic EM algorithm and switch procedure. The comparison of the estimation accuracy between the dynamic EM algorithm with switch procedure and non-dynamic EM algorithm with switch procedure is shown in Table 6. We also performed the simulation with a non-dynamic non-switch procedure for comparison in the Supplementary Materials.

Remark 3.1**.**

The dynamic switch EM algorithm was started with initial ROI $[-10,10]$ cm $\times[-10,10]$ cm $\times[0,10]$ cm and mesh grids $K_{i}=10$, $i=1,2,3$. The ROI was shrunk and the mesh grids were increased by $1$ during iterations. The non-dynamic switch EM algorithm was implemented with the initial ROI and mesh grids throughout all the iterations.

4 Real Data Application $1$

The first real data analysis reports the source localization for the Brain-Controlled Interfaces (BCI) data collected at the Center for Advanced Brian Magnetic Source Imaging (CABMSI) at the Presbyterian University Hospital in Pittsburgh. The data consists of MEG scans of 102 magnetometers recorded at $37000$ milliseconds (ms). During the experiment, the subjects performed a two dimensional center-out task using wrist movement. In the imagined movement task, the subjects were first asked to imagine that they were performing the center-out movement using their wrist. During the overt movement task, subjects controlled a 2-D cursor using their wrist to perform the center-out task. Each trial started after the subject held the cursor in the center for a holding period, followed by a target onset. In order for a trial to be considered successful, the subject needed to move the cursor to the target and hold it there for the duration of the holding period.

The goal of our analysis is to investigate the dynamics of the possible existing sources in the BCI data. Previous work on this data focused on estimating the distribution of the source when the number of sources is assumed to be known [23]. In this section, we consider the data after movement onset from all the magnetometers, and mainly focus on the time varying characteristics of the source location with a dynamically estimated number of sources. We adopt the view-point developed in [24] about the changing number of sources that might exist in the BCI data and exploit the distribution of sources using the estimated number of sources at different stages of the experiment. Throughout the analysis, we assume a unit moment for all possible sources for simplicity. A single sphere head model, with its center $(1.07,0.74,1.65)$ and radius $10.5$, measured in centimeters (cm), was constructed based on the magnetometer positions and head shape information from the BCI data.

Specifically, we have applied the proposed discrete approach to the BCI data through two sub-analyses (short time frame and long time frame): 1a) With the estimated number of sources introduced in [24], we investigate the source distribution in space and time within two selected time windows after the movement onset; 1b) As a previous study has suggested that there are still some active sources present after the movement, we contrast the behavior of the sources for the same selected time windows with no noise estimation; 2) The source estimation for a longer time frame is also reported.

4.1 Activity for Short Time Frame

In 1a, the sources with the estimated number obtained by the Fourier transform (the blue line in Figure 5) were investigated for two selected time windows $[12000,12099]$ ms and $[20000,20099]$ ms, where time window $[12000,12099]$ ms was analyzed for the time varying characteristics of two sources and the time window $[20000,20099]$ ms was selected for the analysis with three sources. In 1b, we considered the same time windows for the source investigation with the estimated number of sources without noise estimation (the red line in Figure 5), which suggested only one and two sources respectively.

• 2.

The Matlab toolbox ‘fieldtrip’ was used to obtain an MNE by searching the entire head. We set the area around the MNE to be the 3D ROI.

• 3.

To implement the proposed discrete approach, we chose mesh grids $K_{i}$, $i=1,2,3$, for $x,y,z$ dimensions of the ROI, respectively. Then, the ROI was subsequently discretized into $K_{1}\cdot K_{2}\cdot K_{3}$ voxels.

• 4.

Motivated by the MNE, we manually set the parameters \mbox{\boldmath\mu}_{0}, \mbox{\boldmath\Sigma}_{0} and $\Sigma$ in the source model, and therefore the movement of the sources was constrained within the ROI all the time. The noise estimation of the selected data was obtained using the Fourier transform.

• 5.

We used the estimates $\hat{A}$ and $\hat{{\bf b}}$ from the converged EM algorithm to calculate the discrete posterior probability distribution of the location parameter $(x,y,z)$ at each time point.

To illustrate the results of the analysis, we now explain the sources distribution (single-source case in 1b and two-source case in 1a) for the time window $[12000,12099]$ ms. For the single-source case, the trajectories of posterior means at each time point during EM iterations are shown in Figure 6 (a), and the posterior means of the converged source are highlighted in Figure 6 (b) at six selected time points. To visualize the time varying characteristics, the target posterior distribution for the location parameters with non-zero probabilities in the 3D ROI are provided. Figure 7 shows the dynamics of the target posterior distribution at the same selected time points. As shown in Figure 7, the posterior means at the selected time points are highlighted using green stars, which vary from $(6,1.73,-1.47)$ cm at the first time point $t=12000$ ms to $(5.66,0.53,-1.81)$ cm at the last time point $t=12099$ ms. Figure 8 represents the marginal posterior distribution at three selected time points, from which we can see the dynamics for the location parameter in each dimension. The dynamics of the two-source case for the same time window have been summarized in Figure 9 to Figure 11 with similar interpretation. The results (two-source case in 1b and three-source case in 1a) for the time window $[20000,20099]$ ms can be found in the Supplementary Materials.

4.2 Activity for Long Time Frame

To investigate the effect of window length in source distribution for the BCI data, we carried out the same analysis for a sequence of two time windows starting at the same time point that was previously studied, but differing in length. In this case, each time window was selected with four different lengths ($200$, $300$, $400$, and $500$ ms, respectively). The distributions of the two sources estimated by the Fourier transform were investigated for four selected time windows $[12000,12k^{*}99]$ ms, where $k^{*}=1,2,3,4$. We chose the ROI and discretized it, similar to the manner that it was done in the previous section, but the noise estimation of the selected data was obtained by the Fourier transform using the measurements within the first $100$ ms. Although the data covariance matrix in noise estimation changes over the length of the window, we have observed that it is necessary to use early observations for stable noise estimation. The posterior distributions for the two sources were obtained for location parameter $(x,y,z)$ at each time point in the selected time windows $[12000,12k^{*}99]$ ms, $k^{*}=1,2,3,4$. To illustrate the results from the four selected time windows, we compared the source distributions with the result from the previous analysis from the time window $[12000,12099]$ ms. Since we considered five time windows with different lengths, we denote $\mathcal{Y}_{[12000,12k^{*}99]}$ as the measurements in the time window $[12000,12k^{*}99]$ ms and \hat{\mbox{\boldmath\Theta}}_{[12000,12k^{*}99]} as the estimate obtained from the measurements $\mathcal{Y}_{[12000,12k^{*}99]}$, for $k^{*}=0,1,\ldots,4$. For $12000\leq t\leq 12099$, we focused on the posterior distributions \mbox{P}({\bf J}_{t}^{\text{\tiny p}}|\mathcal{Y}_{[12000,12k^{*}99]},\hat{\mbox{\boldmath\Theta}}_{[12000,12k^{*}99]}) for $k^{*}=0,1,\ldots,4$. We found that the estimated posterior means of the two sources changed as we increased $k^{*}$ from [math] to $4$. Figure 12 exhibits such changes in source $1$ for only selected time points. The corresponding plot for source $2$ can be found in the Supplementary Materials. On one hand, this result is expected because the selection on time windows should have affected the parameter estimation in EM and the smoothing procedure of the posterior calculation; on the other hand, this serves as an example that the estimation of the source distribution based on more data could result in a significant difference from the short frame. Moreover, we also noted that the changes in the posterior means were subtle in circumstances where $k^{*}$ was increased from [math] to $3$, as shown in the Supplementary Materials. However, when $k^{*}$ was greater than $3$, the marginal posterior means was largely different from those with $k^{*}$ less than 3. This difference seems be caused by a temporary effect. This result matches well with our observation that the filtering result \mbox{P}({\bf J}_{t}^{\text{\tiny p}}|\mathcal{Y}_{[12000,12099]},\hat{\mbox{\boldmath\Theta}}_{[12000,12099]}) might not be a good estimate of the distribution \mbox{P}({\bf J}_{t}^{\text{\tiny p}}|\mathcal{Y}_{[12000,12k^{*}99]},\hat{\mbox{\boldmath\Theta}}_{[12000,12k^{*}99]}) when $k^{*}$ is bigger than $3$. A similar phenomenon was observed for the time window $[20000,20k^{*}99]$ at $k^{*}=3$, where $k^{*}$ was increased from [math] to $4$ with no noise estimation. The details can be found in the Supplementary Materials. We did not find the temporary effect for one-source case in $[12000,12k^{*}99]$ ms and three-source case in $[20000,20k^{*}99]$ ms, for all $k^{*}=0,1,\ldots,4$. All other analyses are in the Supplementary Materials.

5 Real Data Application $2$

For the second real data application, we extend the proposed methodology to a set of EEG recordings under spatial working memory (SWM) task. The SWM reflected in brain activities is often related to relevant brain networks. Thus, source localization studies using the EEG recordings should provide insight into how the brain temporally and spatially response to different SWM loads.

The EEG recordings are collected from a participant under three memory load conditions, which consist of $26$, $27$ and $29$ trials respectively. For each trial, it contains three phases: encoding, retention and probing. During the encoding phase, it began with a cross in the center of the screen. Depending on the load condition, $1$, $3$ or $5$ white dots was/were presented sequentially on the screen. During the retention phase, a fixation cross was then displayed, followed by a red dot presented for the probing phase. The participant was required to indicate whether the red dot appeared at a previously occupied location. The EEG data are processed and down sampled to $250$ Hz. For each trial of EEG recording for the retention phase, it consists of a baseline duration, time before event onset, lasting for $200$ ms and an event duration lasting for the following $3992$ ms.

To investigate the association between EEG responses and brain network in the SWM task, we identify the source distribution at each time step during the event onset period for the retention phase. Since the SWM task often observe sustained negative activity during retention [12], three ROIs associated with fMRI-based deactivation pattern in higher capability group are used for source searching in our study. For each memory load condition, the source distribution is obtained with the estimated number of sources trial by trial. In Figure 13, we plot the source location with significant posterior probability in the three highlighted ROIs (red areas) within the human cortex (grey dots). To distinct brain network under different load conditions, the source distribution from four trials (with $2$-$5$ estimated sources) for each load condition are illustrated column by column. With two or three estimated sources (the first two rows in Figure 13), the source distributions under load $1$ and load $5$ share some similarity that the two/three sources are concentrated in the rightmost ROI, while the source distribution under load $3$ spread out to the other two ROIs. When more sources are presented, the distribution pattern changes. In the last row of Figure 13, it shows the source distribution with five estimated sources. The patterns of the source distributions under load $3$ and load $5$ are similar, with most of the sources distributed in the rightmost ROI. However, none of the sources distributed in the rightmost ROI under load $1$ at the selected time point.

6 Conclusion

The quantification of the source current in a time-varying source model for the MEG data is still a practically urgent problem. Due to the non-linearity of the measurement with respect to the source location parameter, effective methods for dealing with this problem are lacking. Common regularization-based methods mainly focus on estimating the source moment, while other time-varying methods attempt to find estimates of the source location using measurements up to the time point of interest. The latter is encouraged for improvement as using the entire set of measurements is more justifiable for estimating the time-varying source. Meanwhile, these methods are often restricted to a pre-determined number of sources. In MEG, a reliable estimation of the number of sources present in the data, or perhaps at different stages of the experiment, is quite crucial for the whole problem.

With the goal of proposing a framework that allows the flexible estimation of the evolution of the source in MEG, we introduced a discrete approach for calculating the posterior distribution of the source. We emphasized the importance of directly calculating the posterior distribution of the source through a discrete model, and this method differs in the fact that it only samples the intractable continuous target distribution. In this respect, this discrete approach has improved the posterior distribution by providing the probability of the possible sources present in the brain rather than using some approximated samples [23]. In both the single-source case and mutiple-source case, the proposed approach was seen to be more reliable in estimating the source distribution.

We show the sub-optimality of the EM algorithm in estimating the model parameters within the discrete source model. The performance of our proposed approach was examined in some simulated examples, with varying model setups. For the single-source case, the combined procedure seemed to well capture the true time evolution of the sources. For the multi-source model, we have proposed a dynamic procedure and a switch procedure for estimation accuracy and computational efficiency. The proposed approach gave rather satisfactory and consistent results.

In our analysis, we adopted the results of the estimated number of sources in the real MEG/EEG data application from [24], and implemented the source estimation dynamically with and without noise estimation. As suggested in [24], we did not estimate the number of sources for the BCI data from the entire set of data at once, but rather attempted to estimate it from some selected windows. We found the existence of a temporal effect at several selected time points in this data. In fact, we noted that the source distribution estimated from the proposed method and the sequential sampling method in [23] did not differ much in the beginning but started to diverge as more measurements were included. In addition, this phenomenon was only found to be significant for two-source cases in the BCI data. However, the investigation for the EEG data in the SWM task was conducted with the entire measurements after event onset in retention.

To summarize, a reliable estimation of the source distribution depends on a reliable estimation of the number of the sources and the source localization algorithm. With our method, we explore the use of combined approaches in a more advanced form to further examine the evolution of sources.

{supplement}\sname

Supplementary Materials \stitleQuantifying Time-Varying Sources in Magnetoencephalography – A Discrete Approach \slink[doi]COMPLETED BY THE TYPESETTER \sdatatype.pdf \sdescription We include all materials omitted from the main text.

Acknowledgements

The authors thank Professor Rob Kass and his collaborators for sharing their BCI data. The authors would also like to thank collaborator Professor Helen Zhou, at the Duke-NUS medical school, for the spatial working memory (SWM) experiments.

Appendix A

.

Proof of Theorem 2.1. In the $j$-th iteration of the EM algorithm, we are given the estimate \mbox{\boldmath\Theta}^{(j-1)} from the previous iteration. From (2.15), we have

[TABLE]

The update \mbox{\boldmath\Theta}^{(j)} is obtained by maximizing the function L(\mbox{P}(\mathcal{J}^{\text{p}}_{T}|\mathcal{Y}_{T},\mbox{\boldmath\Theta}^{(j-1)}), \mbox{\boldmath\Theta},\mathcal{Y}_{T}), thus we have

[TABLE]

where the second inequality comes from Jensen’s inequality (2.14). Form (A.1) and (A.2), we have that the EM sequence \{\mbox{\boldmath\Theta}^{(j)}\} does not cause a decrease in the log likelihood function. Under the regularity conditions, we refer the poof of Theorem 3 in [21] for the convergence of the EM sequence.

∎

Appendix B

.

Proof of Theorem 2.2. We first consider the case with two sources. From the EM algorithm with the switch procedure and the non-switch procedure, there exists $\delta>0$, such that \|\hat{\mbox{\boldmath\Theta}}_{\text{s}}-\hat{\mbox{\boldmath\Theta}}_{\text{ns}}\|\leq\delta.

In the non-switch procedure, we calculate the posterior distribution of two sources \mbox{P}({\bf J}_{t,1}^{\text{\tiny p}}\in V_{k_{1}},{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}}|\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{ns}}), and can further calculate the marginal posterior distribution

[TABLE]

where $n=1,2$ and $n^{\prime}\neq n$. In the switch procedure, we focus on calculating the marginal posterior distribution \mbox{P}({\bf J}_{t,n}^{\text{\tiny p}}\in V_{k_{n}}|{\bf J}_{t,n^{\prime}}^{\text{\tiny p}}\in V_{k_{n^{\prime}}},\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{s}}), for $n=1,2$, $n^{\prime}\neq n$. For the first source current ${\bf J}^{\text{\tiny p}}_{t,1}$, we aim to calculate

\displaystyle\bigg{|}\mbox{P}(v_{tk_{1}}=1|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{s}})-\sum_{k_{2}=1}^{K}\mbox{P}(v_{tk_{1}}=1,v_{tk_{2}}=1|\mathcal{Y}_{T},\hat{\mbox{\boldmath\Theta}}_{\text{ns}})\bigg{|}

$\displaystyle=$ \displaystyle\bigg{|}\frac{\mbox{P}(v_{tk_{1}}=1,\mathcal{Y}_{t}|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}})\mbox{P}(\mathcal{Y}_{T\setminus t}|v_{tk_{1}}=1,{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}})}{p(\mathcal{Y}_{t}|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}})p(\mathcal{Y}_{T\setminus t}|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}})}-

\displaystyle\quad\quad\sum_{k_{2}=1}^{K}\frac{\mbox{P}(v_{tk_{1}}=1,v_{tk_{2}}=1,\mathcal{Y}_{t}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}})\mbox{P}(\mathcal{Y}_{T\setminus t}|v_{tk_{1}}=1,v_{tk_{2}}=1,\hat{\mbox{\boldmath\Theta}}_{\text{ns}})}{p(\mathcal{Y}_{t}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}})p(\mathcal{Y}_{T\setminus t}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}})}\bigg{|}

$\displaystyle=$ \displaystyle\bigg{|}\alpha^{\text{s}}_{tk_{1}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}})\beta^{\text{s}}_{tk_{1}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}})-\sum_{k_{2}=1}^{K}\alpha^{\text{ns}}_{t,k_{1},k_{2}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}})\beta^{\text{ns}}_{t,k_{1},k_{2}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}})\bigg{|}

(B.1)

,

where we let \alpha^{\text{s}}_{tk_{1}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}}):=\mbox{P}(v_{tk_{1}}=1,\mathcal{Y}_{t}|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}})/p(\mathcal{Y}_{t}|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}}), and \beta^{\text{s}}_{tk_{1}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}}):=\mbox{P}(\mathcal{Y}_{T\setminus t}|v_{tk_{1}}=1,{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}})/p(\mathcal{Y}_{T\setminus t}|{\bf J}_{t,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}}) in the switch procedure, and \alpha^{\text{ns}}_{t,k_{1},k_{2}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}}):=\mbox{P}(v_{tk_{1}}=1,v_{tk_{2}}=1,\mathcal{Y}_{t}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}})/p(\mathcal{Y}_{t}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}}) and \beta^{\text{ns}}_{t,k_{1},k_{2}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}}):=\mbox{P}(\mathcal{Y}_{T\setminus t}|v_{tk_{1}}=1,v_{tk_{2}}=1,\hat{\mbox{\boldmath\Theta}}_{\text{ns}})/p(\mathcal{Y}_{T\setminus t}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}}) in the non-switch procedure. In order to bound (B.1), we analyze $\alpha$’s and $\beta$’s from two procedures recursively. We suppose the same ROI and mesh grids for the discretization $\{V_{k}\}_{k=1}^{K}$.

First, we analyze $\alpha$’s in the forward recursion. For $t=1$, we have

[TABLE]

where c_{1}^{\text{s}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}})=p({\bf Y}_{1}|{\bf J}_{1,2}^{\text{\tiny p}}\in V_{k_{2}},\hat{\mbox{\boldmath\Theta}}_{\text{s}}) in the switch procedure for $1\leq k_{1}\leq K$, and

[TABLE]

where c_{1}^{\text{ns}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}})=p({\bf Y}_{1}|\hat{\mbox{\boldmath\Theta}}_{\text{ns}}) in the non-switch procedure for $1\leq k_{1},k_{2}\leq K$. From regularity condition (C2), there exists positive constants $c$ and $C$, such that c\leq c_{1}^{\text{s}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}}),c_{1}^{\text{ns}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}})\leq C, and we have

[TABLE]

where $c$ is a positive constant. Throughout the proof, we use $c$ as a genetic symbol for positive constant. Since \|\hat{\mbox{\boldmath\Theta}}_{\text{s}}-\hat{\mbox{\boldmath\Theta}}_{\text{ns}}\|\leq\delta, there exist $\varepsilon>0$, such that

[TABLE]

given that \bigg{\|}{\bf c}_{k_{2}^{\ast}}-\sum_{k_{2}=1}^{K}\sum_{k_{1}=1}^{K}\mbox{P}(v_{1k_{1}}=1,v_{1k_{2}}=1|\hat{\mbox{\boldmath\Theta}}_{\text{ns}}){\bf c}_{k_{2}}\bigg{\|}\leq\epsilon, where ${\bf c}_{k_{2}^{\ast}}$ is the location of the second source assumed in the switch procedure. From (B.2) and (B.3), we have

[TABLE]

From item $3$ in Table 1, we have

[TABLE]

and

[TABLE]

for $t=2,\ldots,T$. Therefore, we have

[TABLE]

where (a) follows from (B.2), (b) follows from (B.3) and (B.4), $2\leq t\leq T$, and $1\leq k_{1}\leq K$. Second, we analyze $\beta$’s in the backward recursion. For $t=T$, we initialize \beta^{\text{s}}_{Tk_{1}}(\hat{\mbox{\boldmath\Theta}}_{\text{s}})=\beta^{\text{ns}}_{T,k_{1},k_{2}}(\hat{\mbox{\boldmath\Theta}}_{\text{ns}})=1, for $1\leq k_{1},k_{2}\leq K$. Thus,

[TABLE]

From item $3$ in Table 2, we have

[TABLE]

and

[TABLE]

for $t=T-1,\ldots,1$. Using the same derivation in the forward recursion of (B.5), we have

[TABLE]

Thus, we have

[TABLE]

where $(a)$ follows from (B.1), $(b)$ follows form (B.4), (B.5), (B.6) and (B.7), for $1\leq t\leq T$, $1\leq k_{1}\leq K$. We can obtain the same result for the second source ${\bf J}_{t,2}^{\text{\tiny p}}$. Without loss of generality, we can extend the result to the case with $N$ sources, where $N>2$. ∎