Group level MEG/EEG source imaging via optimal transport: minimum Wasserstein estimates

TL;DR

This paper introduces Minimum Wasserstein Estimates (MWE), a novel group-level MEG/EEG source imaging method that leverages optimal transport to improve localization accuracy by considering cortical geometry and spatial proximity across subjects.

Contribution

The paper proposes a new joint source localization method using optimal transport metrics, enhancing group analysis by accounting for cortical geometry and spatial similarities.

Findings

MWE outperforms standard techniques in simulated data.

MWE accurately localizes neural sources in real MEG data.

Optimal transport-based regularization improves group source estimation.

Abstract

Magnetoencephalography (MEG) and electroencephalogra-phy (EEG) are non-invasive modalities that measure the weak electromagnetic fields generated by neural activity. Inferring the location of the current sources that generated these magnetic fields is an ill-posed inverse problem known as source imaging. When considering a group study, a baseline approach consists in carrying out the estimation of these sources independently for each subject. The ill-posedness of each problem is typically addressed using sparsity promoting regularizations. A straightforward way to define a common pattern for these sources is then to average them. A more advanced alternative relies on a joint localization of sources for all subjects taken together, by enforcing some similarity across all estimated sources. An important advantage of this approach is that it consists in a single estimation in which all measurements are pooled together, making the inverse problem better posed. Such a joint estimation poses however a few challenges, notably the selection of a valid regularizer that can quantify such spatial similarities. We propose in this work a new procedure that can do so while taking into account the geometrical structure of the cortex. We call this procedure Minimum Wasserstein Estimates (MWE). The benefits of this model are twofold. First, joint inference allows to pool together the data of different brain geometries, accumulating more spatial information. Second, MWE are defined through Optimal Transport (OT) metrics which provide a tool to model spatial proximity between cortical sources of different subjects, hence not enforcing identical source location in the group. These benefits allow MWE to be more accurate than standard MEG source localization techniques. To support these claims, we perform source localization on realistic MEG simulations based on forward operators derived from MRI scans. On a visual task dataset, we demonstrate how MWE infer neural patterns similar to functional Magnetic Resonance Imaging (fMRI) maps.

Figures

Equations

$$\mathbf{Y}=\mathbf{B}+\epsilon =\mathbf{Lx}+\epsilon ,$$

$${\mathbf{x}}^{\star }={\mathrm{argmin}}_{\mathbf{x}\in {\mathbb{R}}^p}\Vert \mathbf{Y}-\mathbf{Lx}{\Vert }_2^2.$$

$${\mathbf{x}}^{\star }={\mathrm{argmin}}_{\mathbf{x}\in {\mathbb{R}}^p}\frac{1}{2n}\Vert \mathbf{Y}-\mathbf{Lx}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1,$$

$$\underset{{\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(S)}\in {\mathbb{R}}^p}{\min }\frac{1}{2n}\sum _{s=1}^S\Vert {\mathbf{Y}}^{(s)}-{\mathbf{L}}^{(s)}{\mathbf{x}}^{(s)}{\Vert }_2^2+\Omega ({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(S)}).$$

$$\text{WK}(\mathbf{a},\mathbf{b})=\underset{\begin{array}{c}\mathbf{P}\in {{\mathbb{R}}_{+}}^{p\times p}\\ \mathbf{P}\text{mathds}1=\mathbf{a},{\mathbf{P}}^{\top }\text{mathds}1=\mathbf{b}\end{array}}{\min }⟨\mathbf{P},\mathbf{M}⟩.$$

$$\underset{\mathbf{P}\in {{\mathbb{R}}_{+}}^{p\times p}}{\min }⟨\mathbf{P},\mathbf{M}⟩+\gamma \text{KL}(\mathbf{P}\text{mathds}1|\mathbf{a})+\gamma \text{KL}({\mathbf{P}}^{\top }\text{mathds}1|\mathbf{b}),$$

$$W(\mathbf{a},\mathbf{b})=\underset{\mathbf{P}\in {{\mathbb{R}}_{+}}^{p\times p}}{\min }\epsilon \text{KL}(\mathbf{P}|e^{-\frac{\mathbf{M}}{\epsilon }})+\gamma \text{KL}(\mathbf{P}\text{mathds}1|\mathbf{a})+\gamma \text{KL}({\mathbf{P}}^{\top }\text{mathds}1|\mathbf{b}),$$

$$\tilde{W}(\mathbf{a},\mathbf{b})\stackrel{\mathrm{def}}{=}W({\mathbf{a}}^{+},{\mathbf{b}}^{+})+W({\mathbf{a}}^{-},{\mathbf{b}}^{-}).$$

$${\Omega }_{\text{MTW}}({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(S)})\stackrel{\mathrm{def}}{=}\mu \underset{\bar{\mathbf{x}}\in {\mathbb{R}}^p}{\min }\frac{1}{S}\sum _{s=1}^S\tilde{W}({\mathbf{x}}^{(s)},\bar{\mathbf{x}})+\lambda \Vert {\mathbf{x}}^{(s)}{\Vert }_1,$$

$$\underset{\begin{array}{c}{\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(S)}\in {\mathbb{R}}^p\\ {\sigma }^{(1)},\dots ,{\sigma }^{(S)}\in [{\sigma }_0,+\infty ]\end{array}}{\min }\sum _{s=1}^S\frac{1}{2n{\sigma }^{(s)}}\Vert {\mathbf{Y}}^{(s)}-{\mathbf{L}}^{(s)}{\mathbf{x}}^{(s)}{\Vert }_2^2+\frac{{\sigma }^{(s)}}{2}+{\Omega }_{\text{MTW}}({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(S)}),$$

$${\sigma }^{(s)}\leftarrow \frac{\Vert {\mathbf{Y}}^{(s)}-{\mathbf{L}}^{(s)}{\mathbf{x}}^{(s)}{\Vert }_2}{\sqrt{n}}\wedge {\sigma }_0,$$

$$\underset{\mathbf{x}\in {\mathbb{R}}_{+}^p}{\min }\frac{1}{2n}\Vert \mathbf{Y}-\mathbf{Lx}{\Vert }_2^2+\frac{\mu \gamma }{S}(⟨\mathbf{x},\text{mathds}1⟩-⟨\log (\mathbf{x}),\mathbf{m}⟩)+\lambda \sigma \Vert \mathbf{x}{\Vert }_1,$$

Full text

1 Introduction

Magnetoencephalography (MEG) measures the components of the magnetic field surrounding the head, while Electroencephalography (EEG) measures the electric potential at the surface of the scalp. Both can do so with a temporal resolution of less than a millisecond. Localizing the underlying neural activity on a high resolution grid of the cortex, a problem known as source imaging, is inherently an “ill-posed” linear inverse problem: Indeed, the number of potential sources is larger than the number of MEG and EEG sensors, which implies that, even in the absence of noise, different neural activity patterns could result in the same electromagnetic field measurements.

To limit the set of possible solutions, prior hypotheses on the nature of the source distributions are necessary. The minimum-norm estimates (MNE) for instance are based on $\ell_{2}$ Tikhonov regularization which leads to a linear solution [11]. An $\ell_{1}$ norm penalty was also proposed by [34], modeling the underlying neural pattern as a sparse collection of focal dipolar sources, hence their name “Minimum Current Estimates” (MCE). These methods have inspired a series of contributions in source localization techniques relying on noise normalization [6, 28] to correct for the depth bias [1] or block-sparse norms [30, 10] to leverage the spatio-temporal dynamics of MEG signals. While such techniques have had some success, source estimation in the presence of complex multi-dipole configurations remains a challenge. In this work we aim to leverage the anatomical and functional diversity of multi-subject datasets to improve localization results.

Related work.

This idea of using multi-subject information to improve statistical estimation has been proposed before in the neuroimaging literature. In [20] it is showed that different anatomies across subjects allow for point spread functions that agree on a main activation source but differ elsewhere. Averaging across subjects thereby increases the accuracy of source localization. On fMRI data, [35] proposed a probabilistic dictionary learning model to infer activation maps jointly across a cohort of subjects. A similar idea led [19] to introduce a Bayesian framework to account for functional intersubject variability. To our knowledge, the only contribution formulating the problem as a multi-task regression model employs a Group Lasso with an $\ell_{21}$ block sparse norm [21]. Yet this forces every potential neural source to be either active for all subjects or for none of them.

Contribution.

The assumption of identical functional activity across subjects is clearly not realistic. Here we investigate several multi-task regression models that relax this assumption. One of them is the multi-task Wasserstein (MTW) model [15]. MTW is defined through an Unbalanced Optimal Transport (UOT) metric that promotes support proximity across regression coefficients. However, applying MTW to group level data assumes that the signal-to-noise ratio is the same for all subjects. We propose to build upon MTW and alleviate this problem by inferring estimates of both sources and noise variance for each subject. To do so, we follow similar ideas that lead to the concomitant Lasso [27, 31, 25] or the multi-task Lasso [24].

This paper is organized as follows. Section 2 introduces the multi-task regression source imaging problem. Section 3 presents some background on UOT metrics and explains how MWE are carried out. Section 4 presents the results of our experiments on both simulated and MEG datasets.

Notation.

We denote by $\mathds{1}_{p}$ the vector of ones in $\mathbb{R}^{p}$ and by $\llbracket q\rrbracket\,$ the set $\{1,\ldots,q\}$ for any integer $q\in\mathbb{N}$. The set of vectors in $\mathbb{R}^{p}$ with non-negative (resp. positive) entries is denoted by $\mathbb{R}^{p}_{+}$ (resp. $\mathbb{R}^{p}_{++}$). On matrices, $\log$, $\exp$ and the division operator are applied elementwise. We use $\odot$ for the elementwise multiplication between matrices or vectors. If $\mathbf{X}$ is a matrix, $\mathbf{X}_{i.}$ denotes its $i^{\text{th}}$ row and $\mathbf{X}_{.j}$ its $j^{\text{th}}$ column. We define the Kullback-Leibler (KL) divergence between two positive vectors by $\text{KL}(\mathbf{x},\mathbf{y})=\langle\mathbf{x},\log(\mathbf{x}/\mathbf{y})\rangle+\langle\mathbf{y}-\mathbf{x},\mathds{1}_{p}\rangle$ with the continuous extensions $0\log(0/0)=0$ and $0\log(0)=0$. We also make the convention $\mathbf{x}\neq 0\Rightarrow\text{KL}(\mathbf{x}|0)=+\infty$. The entropy of $\mathbf{x}\in\mathbb{R}^{n}$ is defined as $H(\mathbf{x})=-\langle\mathbf{x},\log(\mathbf{x})-\mathds{1}_{p}\rangle$. The same definition applies for matrices with an element-wise double sum.

2 Source imaging as a multi-task regression problem

We formulate in this section the inverse problem of interest in this paper, and recall how a multi-task formulation can be useful to carry out a joint estimation of all these parameters through regularization.

Source modeling.

Using a volume segmentation of the MRI scan of each subject, the positions of potential sources are constructed as a set of coordinates uniformly distributed on the cortical surface of the gray matter. Moreover, synchronized currents in the apical dendrites of cortical pyramidal neurons are thought to be mostly responsible for MEG signals [26]. Therefore, the dipole orientations are usually constrained to be normal to the cortical surface. We model the current density as a set of focal current dipoles with fixed positions and orientations. The purpose of source localization is to infer their amplitudes. The ensemble of possible candidate dipoles forms the source space.

Forward modeling.

Let $n$ denote the number of sensors (EEG and/or MEG) and $p$ the number of sources. Following Maxwell’s equations, at each time instant, the measurements $\mathbf{B}\in\mathbb{R}^{n}$ are a linear combination of the current density $\mathbf{x}\in\mathbb{R}^{p}:\mathbf{B}=\mathbf{L}\mathbf{x}$. However, we observe noisy measurements $\mathbf{Y}\in\mathbb{R}^{n}$ given by:

[TABLE]

where $\varepsilon$ is the noise vector. The linear forward operator $\mathbf{L}\in\mathbb{R}^{n\times p}$ is called the leadfield or gain matrix, and can be computed by solving Maxwell’s equations using the Boundary element method [12]. Up to a whitening pre-processing step, $\varepsilon$ can be assumed Gaussian distributed $\mathcal{N}(0,\sigma I_{n})$.

Source localization.

Source localization consists in solving in $\mathbf{x}$ the inverse problem (1) which can be cast as a least squares problem:

[TABLE]

Since $n\ll p$, problem (2) is ill-posed and additional constraints on the solution $\mathbf{x}^{\star}$ are necessary. When analyzing evoked responses, one can promote source configurations made of a few focal sources, e.g. using the $\ell_{1}$ norm. This regularization leads to problem (3) called minimum current estimates (MCE), also known in the machine learning community as the Lasso [32].

[TABLE]

where $\lambda>0$ is a tuning hyperparameter.

Common source space.

Here we propose to go beyond the classical pipeline and carry out source localization jointly for $S$ subjects. First, dipole positions (features) must correspond to each other across subjects. To do so, the source space of each subject is mapped to a high resolution average brain using morphing where the sulci and gyri patterns are matched in an auxiliary spherical inflating of each brain surface [8]. The resulting leadfields $\mathbf{L}^{(1)},\dots,\mathbf{L}^{(S)}$ have therefore the same shape $(n\times p)$ with aligned columns.

Multi-task framework.

Jointly estimating the current density $\mathbf{x}^{(s)}$ of each subject $s$ can be expressed as a multi-task regression problem where some coupling prior is assumed on $\mathbf{x}^{(1)},\dots,\mathbf{x}^{(S)}$ through a penalty $\Omega$:

[TABLE]

Following the work of [15], we propose to define $\Omega$ using an UOT metric.

3 Minimum Wassertein Estimates

We start this section with background material on UOT. Consider the finite metric space $(E,d)$ where each element of $E=\{1,\dots,p\}$ corresponds to a vertex of the source space. Let $\mathbf{M}$ be the matrix where $\mathbf{M}_{ij}$ corresponds to the geodesic distance between vertices $i$ and $j$. Kantorovich [16] defined a distance for normalized histograms (probability measures) on $E$. However, it can easily be extended to non-normalized measures by relaxing marginal constraints [4].

Marginal relaxation.

Let $\mathbf{a},\mathbf{b}$ be two normalized histograms on $E$. Assuming that transporting a fraction of mass $\mathbf{P}_{ij}$ from $i$ to $j$ is given by $\mathbf{P}_{ij}\mathbf{M}_{ij}$, the total cost of transport is given by $\langle\mathbf{P},\mathbf{M}\rangle=\sum_{ij}\mathbf{P}_{ij}\mathbf{M}_{ij}$. Minimizing this total cost with respect to $\mathbf{P}$ must be carried out on the set of feasible transport plans with marginals $\mathbf{a}$ and $\mathbf{b}$. The (normalized) Wasserstein-Kantorovich distance reads:

[TABLE]

In practice, if $\mathbf{a}$ and $\mathbf{b}$ are positive and normalized current densities, WK$(\mathbf{a},\mathbf{b})$ will quantify the geodesic distance between their supports along the curved geometry of the cortex. This property makes OT metrics adequate for assessing the proximity of functional patterns across subjects. To allow $\mathbf{a},\mathbf{b}$ to be non-normalized, the marginal constraints in (5) can be relaxed using a KL divergence:

[TABLE]

where $\gamma>0$ is a hyperparameter that enforces a fit to the marginals.

Entropy regularization.

Entropy regularization was introduced by [5] to propose a faster and more robust alternative to the direct resolution of the linear programming problem (5). Formally, this amounts to minimizing the loss $\langle\mathbf{P},\mathbf{M}\rangle-\varepsilon H(\mathbf{P})$ where $\varepsilon>0$ is a tuning hyperparameter. This penalized loss function can be written: $\varepsilon\text{KL}(\mathbf{P},e^{-\frac{\mathbf{M}}{\varepsilon}})$ up to a constant [3]. Combining entropy regularization with marginal relaxation in (6), we get the unbalanced Wasserstein distance as introduced by [4]:

[TABLE]

Generalized Sinkhorn.

Problem (7) can be solved as follows. Let $\mathbf{K}=e^{-\frac{\mathbf{M}}{\varepsilon}}$ and $\psi=\gamma/(\gamma+\epsilon)$. Starting from two vectors $\mathbf{u},\mathbf{v}$ set to $\mathds{1}$ and iterating the scaling operations $\mathbf{u}\leftarrow\left(\mathbf{a}/\mathbf{K}\mathbf{v}\right)^{\psi}$ , $\mathbf{v}\leftarrow\left(\mathbf{b}/\mathbf{K}^{\top}\mathbf{u}\right)^{\psi}$ until convergence, the minimizer of (7) can be computed as $P^{\star}=(\mathbf{u}_{i}\mathbf{K}_{ij}\mathbf{v}_{j})_{i,j\in\llbracket p\rrbracket\,}$. This algorithm is a generalization of the Sinkhorn algorithm [18]. Since it involves matrix-matrix operations, it benefits from parallel hardware, such as GPUs.

Extension to $\mathbb{R}^{p}$.

We extend next the Wasserstein distance to signed measures. We adopt a similar idea to what was suggested in [23, 29, 15] using a decomposition into positive and negative parts, $\mathbf{x}^{(s)}=\mathbf{x}^{(s)_{+}}-\mathbf{x}^{(s)_{-}}$ where $\mathbf{x}^{(s)_{+}}=\max(\mathbf{x}^{(s)_{+}},0)$ and $\mathbf{x}^{(s)_{-}}=\max(-\mathbf{x}^{(s)_{+}},0)$. For any vectors $\mathbf{a},\mathbf{b}\in\mathbb{R}^{p}$, we define the generalized Wasserstein distance as:

[TABLE]

Note that $W(\boldsymbol{0},\boldsymbol{0})=0$ (see [15] for a proof), thus on positive measures $\widetilde{W}=W$. For the sake of convenience, we refer to $\widetilde{W}$ in (8) by the Wasserstein distance, even though it does not verify indiscernability. In practice, this extension allows to compare current dipoles across subjects according to their polarity which could be either towards the deep or superficial layers of the cortex.

The MTW model.

The multi-task Wasserstein model is the specific case of (4) with a penalty $\Omega$ promoting both sparsity and supports’ proximity:

[TABLE]

where $\mu,\lambda\geq 0$ are tuning hyperparameters. The OT term in (9) can be seen as a spatial variance. Indeed, the minimizer $\boldsymbol{\bar{x}}$ corresponds to the Wasserstein barycenter with respect to the distance $\widetilde{W}$.

Minimum Wasserstein Estimates.

One of the drawbacks of MTW is that $\lambda$ is common to all subjects. Indeed, the loss considered in MTW implicitly assumes that the level of noise is the same across subjects. Following the work of [25] on the smoothed concomitant Lasso, we propose to extend MTW by inferring the specific noise standard deviation ${\sigma^{(s)}}$ along with the regression coefficient $\mathbf{x}^{(s)}$ of each subject. This allows to scale the weight of the $\ell_{1}$ according to the level of noise. The Minimum Wasserstein Estimates (MWE) model reads:

[TABLE]

where $\sigma_{0}$ is a pre-defined constant. This lower bound constraint avoids numerical issues when $\lambda\to 0$ and therefore the standard deviation estimate also tends to 0. In practice $\sigma_{0}$ can be set for example using prior knowledge on the variance of the data or as a small fraction of the initial estimate of the standard deviation $\sigma_{0}=\alpha\min_{s}\frac{\|\mathbf{Y}^{(s)}\|}{\sqrt{n}}$. In practice we adopt the second option and set $\alpha=0.01$.

Algorithm.

By combining (7), (8) and (10), we obtain an objective function taking as arguments $\left((\mathbf{x}^{(s)+})_{s},(\mathbf{x}^{(s)-})_{s},(\mathbf{P}^{(s)+})_{s},(\mathbf{P}^{(s)-})_{s},\boldsymbol{\bar{x}}^{+},\boldsymbol{\bar{x}}^{-},(\sigma^{(s)})_{s}\right)$. This function restricted to all parameters except $(\sigma^{(s)})_{s}$ is jointly convex [15]. Moreover, each $\sigma^{(s)}$ is only coupled with the variable $\mathbf{x}^{(s)}$. The restriction on every pair $(\mathbf{x}^{(s)},\sigma^{(s)})$ is also jointly convex [25]. Thus the problem is jointly convex in all its variables. We minimize it by alternating optimization. To justify the convergence of such an algorithm, one needs to notice that the non-smooth $\ell_{1}$ norms in the objective are separable [33]. The update with respect to each $\sigma^{(s)}$ is given by solving the first order optimality condition (Fermat’s rule):

[TABLE]

which also corresponds to the empirical estimator of the standard deviation when the constraint is not active. To update the remaining variables, we follow the same optimization procedure laid out in [15] and adapted to MWE in Algorithm 1. Briefly, let $\mathbf{m}^{(s)+}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{P}^{(s)+}\mathds{1}$ (resp. $\mathbf{m}^{(s)+}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathbf{P}^{(s)+}\mathds{1}$), when minimizing with respect to one $\mathbf{x}^{(s)+}$ (resp. $\mathbf{x}^{(s)-}$), the resulting problem can be written (dropping the exponents for simplicity):

[TABLE]

which can be solved using proximal coordinate descent [7]. Note that the additional inference of a specific $\sigma^{(s)}$ for each subject allows to scale the Lasso penalty depending on their particular level of noise. The final update with respect to $((\mathbf{P}^{(s)+})_{s},(\mathbf{P}^{(s)-})_{s},\boldsymbol{\bar{x}}^{+},\boldsymbol{\bar{x}}^{-})$ can be cast as two Wasserstein barycenter problems, carried out using generalized Sinkhorn iterations [4]. Note that we do not need to compute the transport plans $P^{(s)}$ since inferring every source estimate $\mathbf{x}$ only requires the knowledge of the left marginal $\mathbf{m}=\mathbf{P}\mathds{1}$ which does not require storing $\mathbf{P}$ in memory.

4 Experiments

Benchmarks: Dirty models and Multi-level Lasso.

As discussed in introduction, standard sparse source localization solvers are based on an $\ell_{1}$ norm regularization, applied to the data of each subject independently. We use the independent Lasso estimator as a baseline. We compare MWE to the Group-Lasso estimator [37, 2] which was proposed in this context to promote functional consistency across subjects [21]. It falls in the multi-task framework of (4) where the joint penalty is defined through an $\ell_{21}$ mixed norm $\|\mathbf{X}\|_{21}=\sum_{j=1}^{p}\sqrt{\sum_{s=1}^{S}{\mathbf{x}^{(s)}_{j}}^{2}}$ where $\mathbf{X}\leavevmode\nobreak\ =\leavevmode\nobreak\ (\mathbf{x}^{(s)}_{j})_{(j,s)}\leavevmode\nobreak\ \in\leavevmode\nobreak\ \mathbb{R}^{p\times S}$. We also evaluate the performance of more flexible block sparse models where only a fraction of the source estimates are shared across all tasks: Dirty models [14] and the multivel lasso [22]. In Dirty models source estimates are written as a sum of two parts which are penalized with different norms. One common to all subjects (penalty $\ell_{21}$) and one specific for each subject (penalty $\ell_{1}$). The Multi-level Lasso (MLL) [22] applies the same idea using instead a product decomposition and a Lasso penalty on both parts. We also compare MWE with MTW to evaluate the benefits of inferring noise levels adaptively.

Simulation data and MEG/fMRI datasets.

We use the public dataset DS117 [36] which provides MEG, EEG and fMRI data of 16 healthy subjects to whom were presented images of famous, unfamiliar and scrambled faces. Using the MRI scan of each subject, we compute a source space and its associated leadfield comprising around 2500 sources per hemisphere [9]. Keeping only MEG gradiometer channels, we have $n=204$ observations per subject.

For realistic data simulation, we use the actual leadfields from all subjects, yet restricted to the left hemisphere. We thus have 16 leadfields with $p=2500$. We simulate an inverse solution $\mathbf{x}^{s}$ with $q$ sources ($q$-sparse vector) by randomly selecting one source per label among $q$ pre-defined labels using the aparc.a2009s parcellation of the Destrieux atlas. To model functional consistency, 50% of the subjects share sources at the same locations, the remaining 50% have sources randomly generated in the same labels (see Figure 1). Their amplitudes are taken uniformly between 20 and 30 nAm. Their sign is taken at random with a Bernoulli distribution (0.5) for each label. We simulate $\mathbf{Y}$ using the forward model with a variance matrix $\sigma I_{n}$. We set $\sigma$ so as to have an average signal-to-noise ratio across subjects equal to 4 (SNR$\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\sum_{s=1}^{S}\frac{\|\mathbf{L}^{(s)}\mathbf{x}^{(s)}\|}{S\sigma}$).

We evaluate the performance of all models knowing the ground truth by comparing the best estimates in terms of three metrics: the mean squared error (MSE) to quantify accuracy in amplitude estimation, AUC and a generalized Earth mover distance (EMD) to assess supports estimation. We generalize the PR-AUC (Area under the curve Precision-recall) by defining AUC$(\hat{\mathbf{x}},\mathbf{x}^{\star})=\frac{1}{2}\text{PR-AUC}(\hat{\mathbf{x}}^{+},\mathbf{x}^{\star+})+\frac{1}{2}\text{PR-AUC}(\hat{\mathbf{x}}^{-},\mathbf{x}^{\star-})$ where PR-AUC is computed between the estimated coefficients and the true supports. We compute EMD between normalized values of sources: EMD$(\hat{\mathbf{x}},\mathbf{x}^{\star})=\frac{1}{2}\text{WK}(\frac{\hat{\mathbf{x}}^{+}}{\hat{\mathbf{x}}^{+}\mathds{1}},\frac{\mathbf{x}^{\star+}}{\mathbf{x}^{\star+}\mathds{1}})+\frac{1}{2}\text{WK}(\frac{\hat{\mathbf{x}}^{-}}{\hat{\mathbf{x}}^{-}\mathds{1}},\frac{\mathbf{x}^{\star-}}{\mathbf{x}^{\star-}\mathds{1}})$. Since $\mathbf{M}$ is expressed in centimeters, WK can be seen as an expectation of the geodesic distance between sources. The mean across subjects is reported for all metrics.

Simulation results.

We set the number of sources to 3 and vary the number of subjects under two conditions: (1) using one leadfield for all subjects, (2) using individual leadfields. Each model is fitted on a grid of hyperparameters and the best AUC/MSE/EMD scores are reported. We perform 30 different trials (with different true activations and noise, different common leadfield for condition (1)) and report the mean within a 95% confidence interval in Figure 2.

Various observations can be made. The Group Lasso performs poorly – even compared to independent Lasso – which is expected since sources are not common for all subjects. Non-convexity allows MLL to be very effective with less than 2-4 subjects. Its performance yet degrades with more subjects. OT-based models (MWE and MTW) however benefit from the presence of more subjects by leveraging spatial proximity. They reduce the average error distance from 4 cm (Lasso) to less than 1 cm and reach an AUC of 0.9. One can also observe that the estimation of the noise standard deviation in the MTW model does improve performance. Finally, we can appreciate the improvement of multi-task models when increasing the number subjects, especially when using different leadfield matrices. We argue that the different folding patterns of the cortex across subjects lead to different dipole orientations thereby increasing the chances of source identification.

Results on MEG/fMRI data

The fusiform face area specializes in facial recognition and activates around 170ms after stimulus [17, 13]. To study this response, we perform MEG source localization using Lasso and MWE. We pick the time point with the peak response for each subject within the interval 150-200 ms after visual presentation of famous faces. For both models, we select the smallest $\ell_{1}$ tuning parameter $\lambda$ for which less than 10 sources are active for each subject. Figure 3 shows how UOT regularization favors activation in the ventral pathway of the visual cortex.

The Lasso solutions in Figure 3 show significant differences between subjects. Since no ground truth exists, one could argue that MWE promotes consistency at the expense of individual signatures. To address this concern we compute the standardized fMRI Z-score of the conditions famous vs scrambled faces. We compare Lasso, MWE and fMRI by computing for each subject the ratio largest value in fusiform gyrus / largest absolute value. We report the mean across all subjects in Figure 4. Note that for all subjects, the fMRI Z-score reaches its maximum in the fusiform gyrus, and that MWE regularization leads to more agreement between MEG and fMRI. Figure 5 shows MEG with MWE and fMRI results for subject 2.

Conclusion

We proposed in this work a novel approach to promote functional consistency through a convex model defined using an Unbalanced Optimal Transport regularization. Using a public MEG and fMRI dataset, we presented experiments demonstrating that MWE outperform multi-task sparse models in both amplitude and support estimation. We have shown in these experiments that MWE can close the gap between MEG and fMRI source imaging by gathering data from different subjects.