Prior Variances and Depth Un-Biased Estimators in EEG Focal Source Imaging

TL;DR

This paper introduces a Bayesian method to determine weighting factors in EEG source imaging that eliminate depth bias, leading to more accurate focal source reconstructions.

Contribution

It develops a novel approach to compute weights in sparsity priors for EEG source imaging, ensuring depth-unbiased MAP estimates.

Findings

Proposes Gaussian prior variances for depth-unbiased estimates

Derives approximate weights for sparsity priors based on these variances

Demonstrates improved focal source reconstruction with proposed weights

Abstract

In electroencephalography (EEG) source imaging, the inverse source estimates are depth biased in such a way that their maxima are often close to the sensors. This depth bias can be quantified by inspecting the statistics (mean and co-variance) of these estimates. In this paper, we find weighting factors within a Bayesian framework for the used L1/L2 sparsity prior that the resulting maximum a posterior (MAP) estimates do not favor any particular source location. Due to the lack of an analytical expression for the MAP estimate when this sparsity prior is used, we solve the weights indirectly. First, we calculate the Gaussian prior variances that lead to depth un-biased maximum a posterior (MAP) estimates. Subsequently, we approximate the corresponding weight factors in the sparsity prior based on the solved Gaussian prior variances. Finally, we reconstruct focal source configurations using the sparsity prior with the proposed weights and two other commonly used choices of weights that can be found in literature.

Figures

Tables

	Corresponding author: Alexandra Koulouri
	Institute: University of Münster
	Street: Einsteinstrasse 62
	City: Münster
	Country: Germany
	Email: koulouri@uni-muenster.de

Equations

$$v=Kd+\xi ,$$

$$\pi (d|v)\propto \pi (v|d)\pi (d),$$

$$\pi(v|d)\propto\exp\Big{(}-\frac{1}{2}(v-Kd)^{\mathrm{T}}\Gamma_{\xi}^{-1}(v-Kd)\Big{)}.$$

$$\hat{d}:=\underset{d\in {\mathbb{R}}^{kn}}{\min }\Vert L_{\xi }(Kd-v){\Vert }_2^2-2\ln \pi (d),$$

$$\pi (d)\propto \exp \left(-\frac{1}{2}{d}^{\mathrm{T}}{\Gamma }_d^{-1}d\right)$$

$$\hat{d}=K^{\mathrm{T}}(K{K}^{\mathrm{T}}+{\alpha }^2{\Gamma }_{\xi }{)}^{-1}v.$$

$$\mathbb{E}[\hat{d}]=0\text{mbox}and{\Gamma }_{\hat{d}}=\mathbb{E}[\hat{d}{\hat{d}}^{\mathrm{T}}]=K^{\mathrm{T}}(K{K}^{\mathrm{T}}+{\alpha }^2{\Gamma }_{\xi }{)}^{-1}K.$$

$$\hat{d}={\Gamma }_d{K}^{\mathrm{T}}(K{\Gamma }_d{K}^{\mathrm{T}}+{\alpha }^2{\Gamma }_{\xi }{)}^{-1}v,$$

$${\Gamma }_{\hat{d}}=\mathbb{E}[\hat{d}{\hat{d}}^{\mathrm{T}}]={\Gamma }_d{K}^{\mathrm{T}}(K{\Gamma }_d{K}^{\mathrm{T}}+{\Gamma }_{\xi }{)}^{-1}K{\Gamma }_d.$$

$${\gamma }_d:=\underset{{\gamma }_d}{\min }\Vert \mathrm{diag}({\alpha }^{-2}I-{\Gamma }_{\hat{d}}){\Vert }_2^2.$$

$${\alpha }^2={\gamma }_d^{2(i)}{K}^{(:,i)\mathrm{T}}M{K}^{(:,i)}\text{mbox}fori=1,\dots ,kn,$$

$$\pi (d)\propto \exp \left(-\frac{\alpha }{2}\sum _{i=1}^n{w}_i^r\Vert {\mathrm{d}}_i{\Vert }_2\right)$$

$$\pi (r_i)\propto \exp \left(-\frac{\alpha }{2}{w}_i^r{r}_i\right).$$

$${\gamma }_r^{(i)}=c{\int }_0^{\infty }(r_i-r_{\ast i}{)}^2\exp \left(-\frac{\alpha }{2}{w}_i^r{r}_i\right)d{r}_i=\frac{4}{{\alpha }^2{(w_i^r)}^2}$$

$${\gamma }_r^{(i)}=k{\alpha }^{-2k+2}\left(\prod _{j=1}^k{\gamma }_d^{(i+(j-1)n)}\right){\left(\sum _{j=1}^k{\gamma }_d^{(i+(j-1)n)}\right)}^{-1},$$

$$w_i^r=2\sqrt{\frac{{\alpha }^{2k-4}}{k}\frac{\sum _{j=1}^k{\gamma }_d^{(i+(j-1)n)}}{2\prod _{j=1}^k{\gamma }_d^{(i+(j-1)n)}}}$$

$${\hat{d}}_{\mathrm{MAP}}:=\underset{d}{\min }\Vert v-Kd{\Vert }_2^2+\sum _{i=1}^n\lambda w_i\Vert {\mathrm{d}}_i{\Vert }_{2}w$$

Full text

1 Introduction

In EEG focal source imaging, the goal is to estimate the focal neural activity that arises, for example, during an epileptic seizure using scalp potentials. Based on the distributed source modelling [1], the mapping that connects the dipole moments of $n$ potential source locations to $m$ scalp-potential measurements can be written as

[TABLE]

where $v\in\mathbb{R}^{m}$, $K\in\mathbb{R}^{m\times kn}$ ($m\ll kn$) is the lead field matrix, $k$ is the dimension of the problem (2D or 3D), $d\in\mathbb{R}^{kn}$ is the distributed dipole source configuration and $\xi\sim\mathcal{N}(0,\Gamma_{\xi})$ is the measurement noise.

The ill-posedness of the associated inverse problem requires the use of prior information to obtain stable estimates. One way to solve the problem is to find the estimate of the under-determined linear system that has the minimum norm [2]. However, the minimum norm estimate (MNE) has the property that its maxima can lie only close to the sensors, because the measured scalp potentials can be generated from superficial source configurations with less power than from deep source configurations [3]. Similar source reconstructions can also be obtained with $\ell_{2}$-norm priors. Even if $\ell_{1}$-norm priors are employed the solution consists of several scattered superficial sources [4].

To reduce the depth bias several (often heuristic) approaches have been suggested [5, 6, 3, 7, 8]. The most common approaches are to weight all the sources in the penalty term with the norm of the corresponding column of the lead field matrix [5, 9] or the diagonal elements of the model resolution matrix [10, 11]. Another approach is to use the Bayesian hierarchical modelling [12].

In this paper, our aim is to find, within a Bayesian framework, weights for our sparsity prior such that the resulting posterior estimates do not favor any particular source location or component. Because there is no analytical expression for the MAP estimate when sparsity priors are employed, we propose to solve the weights indirectly. We first quantify the depth bias of the maximum a posterior (MAP) estimates when an i.i.d. Gaussian prior is employed by inspecting the statistics of the MAP estimates. Next, we calculate the Gaussian prior variances that ensure depth un-biased solutions by equalizing the variances in the covariance matrix of the MAP estimates. Finally, we approximate the corresponding weighting factors in the sparsity prior using the solved Gaussian prior variances. We demonstrate the feasibility of our approach by simulating focal brain activity with finite element (FE) simulations. In the reconstructions, we employ the weighted $\ell_{1}/\ell_{2}$ sparsity prior and we compare the results obtained using our proposed weights with the reconstructions based on two other commonly used choices of depth weights.

2 Theory

2.1 Bayesian Inversion

In the Bayesian framework, the inverse solution is the posterior density of the Bayes formula

[TABLE]

where $\pi(v|d)$ is the likelihood and $\pi(d)$ the prior. From Equation (1), the likelihood is

[TABLE]

The MAP estimate of the reconstructions is [13],

[TABLE]

where $L_{\xi}$ comes from the Cholesky factorization of $\Gamma_{\xi}$.

2.2 Gaussian prior

Let us consider a Gaussian prior

[TABLE]

that does not have depth weights i.e. the covariance matrix is $\Gamma_{d}=\alpha^{-2}I$ where $I$ is the identity matrix and $\alpha^{2}$ a scaling parameter. In this case, the MAP estimate is [2]

[TABLE]

From variational point of view, this MAP estimate coincides with Tikhonov regularization and thus yields to a harmonic solution that attains its maximum at the boundary [4]. This can also be explained statistically by analyzing the expectation value and covariance of the MAP estimates. Theses values can be estimated by sampling or by using the analytical expressions

[TABLE]

Figure 1-A shows how the values of the diagonal elements (variances) of $\Gamma_{\hat{d}}$ decrease almost quadratically with respect to depth. The zero expectation values and the very low variances associated with the deep locations imply that the deep sources are very unlikely to be reconstructed. Thus, this MAP estimator is biased with respect to depth and favors sources close to the sensors.

In this paper, our aim is to determine such prior variances that the resulting MAP estimates do not favor any particular source location or component over other i.e. the variances of the MAP estimates are equal.

We start by postulating that this prior covariance matrix is diagonal $\Gamma_{d}=\alpha^{-2}\mathrm{diag}(\gamma_{d}^{(i)})$ for $i=1,\ldots,kn$. The MAP estimate corresponding to this prior is

[TABLE]

and the covariance of the MAP estimates becomes

[TABLE]

In a similar way as in [14], we estimate the prior variances by minimizing

[TABLE]

This results in solving a set of non-linear equations

[TABLE]

where $M=(\Gamma_{{\xi}}+K\Gamma_{{d}}K^{\mathrm{T}})^{-1}$ and ${K}^{(:,i)}$ is the $i^{th}$ column.

Figure 1-B shows that with these prior variances the diagonal elements of $\Gamma_{\hat{d}}$ will be equal, or in other words, the corresponding MAP estimator is depth unbiased. Moreover, Figure 1-C depicts the diagonal elements of the posterior covariance $\Gamma_{d|v}=(K^{\mathrm{T}}\Gamma_{\xi}^{-1}K+\Gamma_{d}^{-1})^{-1}$ obtained based on the estimated prior and Figure 1-D shows two corresponding marginal posterior densities of two different locations. We can observe that the posterior dipole variances increase with respect to depth. Qualitatively, this means that in the estimated source configurations the deep sources are allowed to have higher strengths than the superficial sources, and therefore, the solutions can attain their maximum also deeper in the brain (and not only close to the sensors).

2.3 $\ell_{1}/\ell_{2}$- norm sparsity prior

In this paper, we consider sparse focal source reconstructions and therefore, we employ the $\ell 1/\ell_{2}$-norm prior

[TABLE]

where $d_{i}=(d_{ix},d_{iy},d_{iz})$, $\|d_{i}\|_{2}=\sqrt{d_{ix}^{2}+d_{iy}^{2}+d_{iz}^{2}}$ is the strength of the source at location $i$ and $w_{i}^{r}$ are the weights. For short, we denote the dipole strength at location $i$ as ${r}_{i}=\|\mathrm{d}_{i}\|_{2}$ and

[TABLE]

The variance of $\pi(r_{i})$ is

[TABLE]

where $r_{*i}=c\int_{0}^{\infty}r_{i}\exp{\left(-0.5\alpha w_{i}^{r}r_{i}\right)}~{}dr_{i}=\frac{4c}{\alpha^{2}{(w_{i}^{r})}^{2}}$ and $c=0.5\alpha w_{i}^{r}$ because $\int_{0}^{\infty}c\exp{\left(-0.5\alpha w_{i}^{r}r_{i}\right)}~{}dr_{i}=1$.

We calculate $\gamma_{r}^{(i)}$ at location $i$ with the help of the corresponding Gaussian variances $\gamma_{d}^{(i+(j-1)n)}$ as

[TABLE]

where $j=1,\ldots,k$ and $k$ is the dimension of the problem. This choice ensures that $\gamma_{r}^{(i)}$ is roughly the average of the dipole component variances when the variances of the components are similar and that $\gamma_{r}^{(i)}$ is close to the lowest dipole component variance when the variances have large differences. Finally, from Equation (14) and (15) we calculate the weights

[TABLE]

The estimated Gaussian variances and the corresponding weights of the $\ell_{1}/\ell_{2}$-norm prior are shown in Figure 2.

3 Materials and methods

We study the proposed weights by simulating focal deep sources in the gray matter of a 2D FE head model. The head model consisted of five compartments with conductivities (in S/m) equal to $0.33$ for the scalp, $0.015$ for the skull, $1.76/0.016/0.33$ for the cerebral spinal fluid, gray matter and white matter [15], respectively. The potential measurements $v$ were obtained from 32 point sensors equally spaced around the boundary. For the forward and the inverse computations, we use two meshes with 2342 and 1236 nodes, respectively.

The MAP estimate of the dipole configuration with sparsity constraint is

[TABLE]

where $\lambda$ is a tuning parameter. The minimization is performed by using the interior point method [16] with Bregman iterations [17]. The performance of the proposed weights, $w_{i}^{r}$, from Equation (16), was compared with two other commonly used weights: first, the MNE resolution weights given by $w_{i}^{\mathrm{MNE}}=\sqrt{1/k\sum_{j=1}^{k}R^{(i,i+(j-1)n)}}$, where $R^{(i,i)}=\mathrm{diag}(K^{\mathrm{T}}(KK^{\mathrm{T}}+\Gamma_{\xi})^{-1}K)$ [11] and second, the normalized maximum sensor responses $w^{\mathrm{M}SR}_{i}=g_{i}/\max{(g_{i})}$, where $g_{i}=\max_{l=1:m}\left(\|1/k\sum_{j=1}^{k}K^{(l,i+(j-1)n)}\|_{2}\right)$ [18]. To access the ground truth, we consider measurements with high signal to noise ratio, $\mathrm{SNR}=60$dB. For the quantitative comparison of the results we employ the earth mover’s distance (EMD) [19].

4 Results and discussion

We demonstrate the performance of the different weights using three test cases with one and two dipole sources. In Figure 3, the small images on the left hand side show the true dipoles, the location is marked with blue circles and the orientations with small blue lines. The remaining images, starting from left, show the reconstruction when $w_{i}^{\mathrm{MNE}}$, $w_{i}^{\mathrm{MSR}}$ and $w_{i}^{r}$ are used as weights, respectively. The blue marker $\mathrm{x}$ shows the locations of true sources. The MAP estimates were computed by solving Equation (17).

All the tested weights give feasible reconstructions. However, we note that the proposed weights $w_{i}^{r}$ give the least scattered results and work the best in the single focal source cases. For the two source case, all the weights give roughly similar reconstructions and EMD values.

5 Conclusion and future work

We have demonstrated that the proposed depth weights with the $\ell_{1}/\ell_{2}$ sparsity prior give better reconstruction compared to two commonly used weights when single deep sources are studied. Our proposed approach has the benefit that it does not require using hyper-parameter models that would involve extensive sampling due to the lack of an analytical expression for the MAP estimate when the $\ell_{1}/\ell_{2}$ prior is used. In the future, Monte Carlo simulations will be carried out in a 3D head model to analyze the distribution of the MAP estimates reconstructed by using the $\ell_{1}/\ell_{2}$ sparsity prior with the proposed weights.

Conflict of interest

The authors declare that they have no conflict of interest.