Bayesian multi--dipole localization and uncertainty quantification from simultaneous EEG and MEG recordings

TL;DR

This paper presents a Bayesian method using sequential Monte Carlo to accurately localize multiple brain sources from combined EEG and MEG data, reducing uncertainty compared to single-modality analysis.

Contribution

It introduces a novel Bayesian multi-dipole localization approach that effectively integrates EEG and MEG data for improved source estimation.

Findings

Combined data improves localization of elusive sources.

Posterior distribution variance is lower with combined data.

Method successfully characterizes the number and locations of dipoles.

Abstract

We deal with estimation of multiple dipoles from combined MEG and EEG time--series. We use a sequential Monte Carlo algorithm to characterize the posterior distribution of the number of dipoles and their locations. By considering three test cases, we show that using the combined data the method can localize sources that are not easily (or not at all) visible with either of the two individual data alone. In addition, the posterior distribution from combined data exhibits a lower variance, i.e. lower uncertainty, than the posterior from single device.

Equations

$$y_t^{\text{(E/M)}}=\sum _{n=1}^N{G}^{\text{(E/M)}}(r^n)\cdot q_t^n+n_t^{\text{(E/M)}}$$

$$p(x|y_{1:T}^{\text{(E/M)}})=\frac{p(y_{1:T}^{\text{(E/M)}}|x)p(x)}{p(\text{/}{y}_{1:T}^{\text{(E/M)}}\text{/})}.$$

$$p(y_{1:T}^{\text{(E/M)}}|x)=\prod _t\mathcal{N}\left(y_t-\sum _{n=1}^N{G}^{\text{(E/M)}}(r^n)\cdot q_t^n;0,{\Sigma }^{\text{(E/M)}}\right),$$

$$p(x)=p(N)\prod _{n=1}^N\left(p(r^n)\prod _{t=1}^{T}p(q_t^n)\right)$$

$$p(x|y_{1:T}^{\text{(E/M)}})=p(q_{1:T}^{1:N}|r^{1:N},N,y_{1:T}^{\text{(E/M)}})p(N,r^{1:N}|y_{1:T}^{\text{(E/M)}});$$

$$p_k(N,r^{1:N}|y_{1:T}^{\text{(E/M)}})\propto p(y_{1:T}^{\text{(E/M)}}|N,r^{1:N}{)}^{\alpha (k)}p(N,r^{1:N}),$$

$$\mathbb{P}(N|y_{1:T}^{\text{(E/M)}})=\sum _{i=1}^I{w}_i\delta (N,N_i)$$

$$\mathbb{P}(r|y_{1:T}^{\text{(E/M)}},\hat{N})=\sum _{i=1}^I{w}_i\delta (\hat{N},N_i)\sum _{k=1}^{N_i}\delta (r,r_i^k);$$

$$y_t^{\text{(E+M)}}=\sum _{n=1}^N{G}^{\text{(E+M)}}(r^n)\cdot q_t^n+n_t^{\text{(E+M)}},$$

$$\Gamma^{\text{\hskip 0.8pt\tiny$(\mathbf{E+M})$}}=\left(\begin{array}[]{cc}\Gamma^{\text{\hskip 0.8pt\tiny$(E)$}}&0\\ 0&\Gamma^{\text{\hskip 0.8pt\tiny$(M)$}}\end{array}\right)\ .$$

Full text

1 Introduction

Electroencephalography (EEG) and Magnetoencephalography (MEG) provide non–invasive measurements of brain activity with high temporal resolution. The spatial resolution of these techniques, however, is partially limited by the ill–posedness of the inverse problem. In particular, the two instruments are almost insensitive to currents orientated either radially (MEG) or tangentially (EEG). Therefore EEG and MEG recordings provide complementary information on the underlying neural activity [1, 2].

For this reason, in the last two decades several methods have been proposed that combine EEG and MEG measurements using standard linear inversion [3, 4], sparse linear inversion [5], empirical Bayes with multiple sparse priors [6], maximum entropy on the mean [7], beamformers [8]. However, most of the existing studies that use simultaneous EEG and MEG data rely on a distributed source model, rather than on the dipolar one. Recently we have been working on sequential Monte Carlo methods for approximating the posterior distribution of a set of current dipoles from M/EEG time–series [9, 10, 11, 12]. Here we consider a static multi–dipole model, in which the number of dipoles and their locations do not change with time. This assumption is coherent with the neurophysiological interpretation of a current dipole as the activity of a neural population; in addition, it is commonly used in multi–dipole modeling of M/EEG data. We slightly modify the Sequential Monte Carlo algorithm developed in [11] in order to approximate the posterior distribution of the unknown parameters given simultaneous EEG and MEG recordings and show the advantages of the use of combined MEG and EEG data. Indeed, not only localization of purely radial/tangential sources becomes feasible, but also the uncertainty on the source location gets smaller.

The article is organized as follows. In Section II we review the main features of the algorithm originally proposed in [11] for MEG data, and indicate how it can be modified to use simultaneous MEG and EEG data. In Section III we show three examples of application to synthetic time series. In Section IV we offer a brief discussion and our conclusions.

2 Bayesian multi–dipole modeling from simultaneous EEG and MEG data

2.1 A Bayesian approach to multi–dipole modeling

Let y^{\text{\hskip 0.8pt\tiny(E/M)}}_{1:T}=\left\{y_{t}^{\text{\hskip 0.8pt\tiny(E/M)}}\right\}_{t=1,\dots,T} be the recorded electric potential ($E$) or magnetic field ($M$), sampled at discrete time instants $t$. In the classic multi–dipole model of M/EEG data, it is assumed that the recordings have been produced by a small set of point–like currents, named current dipoles, each of which is defined by its location $r$ and dipole moment $q$. The relationship between the measured data and the dipole parameters is given by

[TABLE]

where: $N$ is the number of dipoles; G^{\text{\hskip 0.8pt\tiny(E/M)}}(r), the leadfield, is a known and non–linear function of the dipole location; n^{\text{\hskip 0.8pt\tiny(E/M)}}_{t} is zero–mean Gaussian noise, with covariance matrix \Sigma^{\text{\hskip 0.8pt\tiny(E/M)}}. We remark that in our model the number of dipoles and their location do not change with time, while dipole moments do.

In a statistical Bayesian approach, one is interested in characterizing the posterior distribution for the unkonwn variables $x=\left\{N,r^{1},q^{1}_{1:T},\dots,r^{N},q^{N}_{1:T}\hskip 0.8pt\right\}$, conditioned on the data:

[TABLE]

In (2), the likelihood function is the following Gaussian

[TABLE]

while the prior distribution is

[TABLE]

where: $p(N)$ is chosen as a Poisson distribution with mean value $\lambda<1$, so as to favour low–cardinality models; $p(r^{\hskip 0.8ptn})$ is uniform in the brain volume; $p(q^{\hskip 0.8ptn}_{t})$ is Gaussian. In the present case, it is useful to split the posterior distribution as follows:

[TABLE]

indeed, under Gaussian assumptions for the prior $p(q^{1:N}_{1:T})$, the first bit at the right–hand side is Gaussian and can be computed analytically, while the second one can be approximated by means of Monte Carlo sampling.

2.2 Adaptive sequential Monte Carlo sampling

Sequential Monte Carlo algorithm. We use a Monte Carlo technique to approximate the posterior distribution for the number and location of dipoles, defined in equation (5). Specifically, we adopt a recently developed class of algorithms known as Sequential Monte Carlo samplers [13] that behave very efficiently in the case of complex posteriors. The two main ideas underlying this class of algorithms are:

• •

instead of trying to sample the posterior directly, they reach the posterior distribution smoothly, through a tempering sequence of auxiliary distributions. Here we adopt the sequence defined as

[TABLE]

with $\alpha(1)=0\hskip 0.8pt$, $\alpha(K)=1$ and $\alpha(1)<\alpha(2)<\ldots<\alpha(K)$, being $K$ the number of iterations;

• •

they combine Importance Sampling (IS) and Markov Chain Monte Carlo (MCMC) [14]: first, they draw a given number $I$ of uniformly weighted samples from the prior distribution; then, they let each sample evolve, at every iteration, according to a Markov Chain Monte Carlo transition kernel (that is $p_{k+1}$–invariant), and re–weight the samples appropriately. Such combination of IS and MCMC crucially combines the benefits of IS (independent samples, global exploration of the state–space) with those of MCMC (correlated samples, good local exploration of the state–space).

At each iteration, the algorithm gets a sample set that is distributed according to p_{k}(N,r^{1:N}\hskip 0.8pt|\hskip 0.8pty_{1:T}^{\text{\hskip 0.8pt\tiny(E/M)}}); sample points can therefore account gradually for the information content of the data. Importantly, the sequence of values of the exponent (and therefore the total number of iterations) needs not to be established a priori; instead, it can be adaptively determined at run–time, based on how well the sample set $\left\{(N_{i},r^{\hskip 0.8pt1:N}_{i}\hskip 0.8pt)\right\}_{i=1,...,I}\hskip 0.8pt$ approximates p_{k}(N,r^{\hskip 0.8pt1:N}\hskip 0.8pt|\hskip 0.8pty_{1:T}^{\text{\hskip 0.8pt\tiny(E/M)}})\hskip 0.8pt.

Estimates. The approximated posterior distribution contains information on multiple alternative models; therefore in order to produce a sensible map, we restrict our attention to the most probable model, by first estimating the posterior probabilities for the number of sources; i.e. we compute the \arg\max\mathbb{P}(N\hskip 0.8pt|\hskip 0.8pty_{1:T}^{\text{\hskip 0.8pt\tiny(E/M)}}), with

[TABLE]

being $w^{i}$ and $N^{i}$ the weight and the number of dipoles of the i–-th sample respectively. Subsequently, denoting by $\hat{N}\hskip 0.8pt$ the most probable number of dipoles, for each voxel $r$, we compute

[TABLE]

this quantity represents the posterior probability of there being a dipole in that voxel; on the other hand, the integral of \mathbb{P}(r\hskip 0.8pt|\hskip 0.8pty_{1:T}^{\text{\hskip 0.8pt\tiny(E/M)}},\hskip 0.8pt\hat{N}) over a volume can be interpreted as the expected number of sources in the integration volume.

2.3 Combining MEG and EEG data

Here we aim at characterizing the posterior distribution for the unknowns given simultaneous recordings of the electric potential and magnetic field. In fact, there will be no formal difference with the single–device scenario. If we concatenate the data and define y_{t}^{\text{\hskip 0.8pt\tiny(E+M)}}=(\/y_{t}^{\text{\hskip 0.8pt\tiny(E)}},\,y_{t}^{\text{\hskip 0.8pt\tiny(M)}}\/), and do the same for the leadfield matrices G^{\text{\hskip 0.8pt\tiny(E+M)}}=[\hskip 0.8ptG^{\text{\hskip 0.8pt\tiny(E)}}\hskip 0.8ptG^{\text{\hskip 0.8pt\tiny(M)}}]\hskip 0.8pt, and the noise vectors n^{\text{\hskip 0.8pt\tiny(E+M)}}_{t}=(\/n^{\text{\hskip 0.8pt\tiny(E)}}_{t},\/n^{\text{\hskip 0.8pt\tiny(M)}}_{t}\/), then we can rewrite equation (1) as:

[TABLE]

where, of course, the underlying hypothesis is that the dipoles producing the measured EEG and MEG data are the same. Then one can go on with exactly the same discussion as before. The key point is to set the noise covariance matrix properly: in the present case, we will be using

[TABLE]

3 Results

We consider three exemplar test cases: two single–dipole cases and one three–dipole case. The leadfield matrices were generated with the BESA Research software. The EEG leadfield contains 32 sensors, the MEG leadfield contains 122 sensors. The source space is discretized with 18,023 voxels (possible dipole locations). The electrical conductivity inside the head was modeled as piecewise constant and homogeneous; a realistic head geometry was used; the leadfield calculation was performed using a boundary element method. Synthetic MEG and EEG data were generated by choosing dipole location(s) and orientation(s), and then assigning dipoles a bell–shaped source time course. White Gaussian noise was added to the synthetic data to produce a Signal–to–Noise Ratio comparable to that of experimental data. The MEG synthetic data generated for the three–dipole case are shown in Figure 1.

In the two single–dipole cases, dipole orientations were selected to make the dipole tangential, in one case, and radial, in the other case, so as to make localization from MEG or EEG only, respectively, particularly challenging.

For each test dataset, we applied the SMC algorithm to the EEG data only, to the MEG data only and to the combined EEG+MEG data. In Figures 2–4 we show the localization results: we highlight the points having empirical posterior probability, as computed by (8) above the threshold $P_{th}=1/1000$. Red, blue and green refer to the results obtained from the EEG data only, from MEG data only, and from combined EEG+MEG, respectively.

As expected, even single dipole localization is challenging, with M/EEG only, when the dipolar sources are tangential/radial. For the tangential source (Fig. 2), the posterior from EEG data is smeared and indicates several possible distinct locations. For the radial source (Fig. 3), the posterior from MEG data concentrates in a deep and central area of the brain, confirming that the signal does not allow proper localization of the source. However, in both cases the posterior from the combined data is able to localize the dipolar source. In order to quantify the improvement, here we compute the average distance from the true dipole location $r^{*}$, i.e.: $\Delta=E_{\mathbb{P}(r|y_{1:T},\hat{N})}\|r-r^{*}\|$. We notice that this quantity contains information on both accuracy of the reconstruction and variance of the distribution. For the tangential dipole, we get $\Delta^{(E)}=28$ mm, $\Delta^{(M)}=8$ mm and $\Delta^{(E+M)}=0.4$ mm; for the radial dipole, we get $\Delta^{(E)}=7$ mm, $\Delta^{(M)}=50$ mm and $\Delta^{(E+M)}=0.4$ mm, again. Remarkably, in both cases the results from the combined measurement are better than those from the best single–measurement case.

For the three–dipole case (Fig. 4), the EEG posterior indicates two sources, none of which corresponds exactly to one of the three true dipoles. The MEG localization is quite accurate, however, also in this case the combined localization is better, featuring a smaller uncertainty. Quantitative evaluation of the performance is more delicate in this case, and is left for a more thorough study.

4 Discussion and Conclusions

The preliminary results contained in this article confirm that the combination of EEG and MEG data can provide better localization with respect to the use of just one device. Interestingly, the use of a Bayesian approach in this context allows to characterize the performance improvement not only in terms of missing sources or localization error, but also in terms of uncertainty of the estimated location, coded in the variance of the posterior probability distribution.

From a statistical perspective, the results of this paper cannot be considered as a fair comparison between EEG, MEG and EEG+MEG, because the number of sensors considered in the three cases is not the same. At the same time, the results can be considered of practical relevance, inasmuch as the comparison has been performed using typical sensor configurations, as available in many research centers and clinics.

Future work will be devoted to investigate this whole problem more systematically, by extending the number of test cases and providing quantitative measures of improvements.

Conflict of Interest

The authors declare that they have no conflict of interest.

Acknowledgements

BESA GmbH, and particularly Andre Waelkens, is kindly acknowledged for providing the leadfield matrix.