# MEM-diffusion MRI framework to solve MEEG inverse problem

Brahim Belaoucha, Jean-Marc Lina, Maureen Clerc, Théodore Papadopoulos

► To cite this version:

Brahim Belaoucha, Jean-Marc Lina, Maureen Clerc, Théodore Papadopoulos. MEM-diffusion MRI framework to solve MEEG inverse problem. 2015 23rd European Signal Processing Conference (EU-SIPCO), Aug 2015, Nice, France. hal-01207165

HAL Id: hal-01207165

https://inria.hal.science/hal-01207165

Submitted on 28 Apr 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# MEM-diffusion MRI framework to solve MEEG inverse problem

Brahim Belaoucha†,‡

Jean-Marc Lina *,**

Maureen Clerc †

Théodore Papadopoulos †

† Athena Project-Team, Inria Sophia Antipolis-Méditerranée, France

‡ Université Côte d'Azur, Nice, France

* Centre de Recherches Mathématique, Montréal, Canada

**Département de Génie Electrique, ETS, Montréal, Canada

## ABSTRACT

In this paper, we present a framework to fuse information coming from diffusion magnetic resonance imaging (dMRI) with Magnetoencephalography (MEG)/Electroencephalography (EEG) measurements to reconstruct the activation on the cortical surface. The MEG/EEG inverse-problem is solved by the Maximum Entropy on the Mean (MEM) principle and by assuming that the sources inside each cortical region follow Normal distribution. These regions are obtained using dMRI and assumed to be functionally independent. The source reconstruction framework presented in this work is tested using synthetic and real data. The activated regions for the real data is consistent with the literature about the face recognition and processing network.

**Index Terms**— MEG, EEG, dMRI, source reconstruction, parcellation, MEM.

## 1. INTRODUCTION

Magneto and electroencephalography (MEEG) are non-invasive modalities which provide an insight on the temporal succession of cognitive processes. Finding the activation is an underdetermined problem due to the large number of unknowns i.e. dipoles (distributed sources) with compared to the number of measurements. A solution is obtained by assuming a prior on the solution. Different prior choices of lead to different approaches [1]. They can be divided into two categories: linear and nonlinear approaches. Linear approaches overestimate the sources' intensities. Nonlinear approaches, like MEM, do not suffer from this limitation.

Physiologically, the brain is organized in functional areas [2]. In the previous MEM source estimation papers, functional imaging (fMRI) and the Multivariate Source Localization (MSP) [3] have been used to parcellate the cortical surface. However, these approaches give different parcelations per subject depending on the functional data. dMRI is a non-invasive imaging modality which gives the diffusion of water in the tissues, as a result it reveals the fiber structure. In this work, we define a unique cortical parcellation for each subject that depends only on the brain anatomy. We use dMRI

information to segment the white/gray (W/G) matter interface into some correlate of functional areas as in Anwander et al. [4], but we extend this approach to work on the whole cortex. We then integrate the cortical surface parcellation in the MEM framework [5] presented in section 2. The sources inside each region are assumed to follow a Normal distribution, and for simplicity we consider them functionally independent. We evaluate the accuracy of the source estimates using synthetic and real MEEG measurements obtained from a face recognition task. The resulting active regions of the MEM-dMRI approach is compared to the Minimum norm estimate (MNE) [1] and what can be found in literature about face processing and recognition network.

## 2. MEM-DMRI FRAMEWORK

### 2.1. MEM

For each task, let $M \in \mathbb{R}^{C \times N \times T}$ ($C$ number of sensors, $N$ number of task repetitions, $T$ the time window length) denote the measurements, (magnetic field for MEG, electrode potential for EEG). $R \in \mathbb{R}^{S \times N \times T}$ ($S$ number of sources) be the source intensities, and $G \in \mathbb{R}^{C \times S}$ be the lead field matrix. $m \in \mathbb{R}^{C \times T}$ is the averaged signal over repetitions of $M$ i.e. $m = E[M]$ and can be written as:

$$m = GE[R] + \varepsilon$$

where $\varepsilon$ is a zero mean Gaussian measurement noise.

The MEM framework maximizes the information coming from the measurements, by maximizing the entropy of the mean signal. Let's assume that the dipole intensities follow the probability law; $p(r) = f(r)\mu(r)$, where $\mu$ is the reference probability distribution of all dipoles (denoted by $r$) in the cortical surface without any MEG/EEG measurements. Let's consider the Kullback-Leibler's pseudo-distance between the reference and the real distribution of the dipoles intensities [5]

$$S_{\mu}(p) = - \int p \log \frac{p}{\mu} \qquad (1)$$

The probability law, $p$, is chosen to maximize the information brought by the measurements. This yields to the definition of
the following objective function L to be minimized:

$$L(p, \lambda, \lambda_0) = -S_{\mu}(p) + \lambda^{t}(m - Gr) + \lambda_0(1 - \int p) \quad (2)$$

where (·) denotes the transpose. The dipole intensity is given by [5]:

$$\tilde{r} = \nabla F_{\mu}(X)|_{X=G^{t}\tilde{\lambda}} \quad (3)$$

where $F_μ$ is the log of the partition function associated to the distribution $μ$

$$F_{\mu}(G^{t}\lambda) = \log \int \exp(\lambda^{t}Gr)\mu(r)$$

and $\tilde{\lambda}$ is the vector defined by;

$$\tilde{\lambda} = \operatorname*{argmax}_{\lambda}(\lambda^{t}m - F_{\mu}) \quad (4)$$

In the presence of noise with the probability measure $d\mu_n$, equation (4) is still valid, we need just to replace $F_μ(G^t\tilde{\lambda})$ by $F_μ(G^t\tilde{\lambda}) + F_{d\mu_n}(\tilde{\lambda})$ [5].

## 2.2. MEM-dMRI

fMRI and MSP were used to parcellate the cortical surface and used in the MEM framework. They give different parcel- lations per subject depending on the functional data. In this paper, we use dMRI to parcellate (see section 3.1) the corti- cal surface into K cortical parcels, {P1, ..., P_K}. The sources in each region are assumed to follow a Normal distribution. Each parcel has two possible states, i.e. for a region i, $S_i \in \{0, 1\}$, which represents inactive and active regions respec- tively. Using this information, the reference probability for independent regions can be found as explained in [5]:

$$\mu(r) = \prod_{k=1}^{K} \sum_{S_k=0}^{1} \pi(S_k) \mu(r_k|S_k) \quad (5)$$

where $\pi$ is the joint probability law of $S_k$ and $r_k$ is the sources inside the $k^{th}$ region. Using this definition and (3), the MEM estimate is [5]:

$$\tilde{r}_{i,k} = \sum_{S} \tilde{\pi}(S) \nabla_{i} F_{S_k}(X)|_{X=G_k^{t} \tilde{\lambda}} \quad (6)$$

where $\nabla_i F_{S_k}(X)|_{X=G_k^t \tilde{\lambda}}$ is the estimation of $r_i$, a dipole i in the parcel $P_k$ for a given state $S_k$, the summation over $S$ runs over all the configuration of $S = \{S_1, S_2, ..., S_k\}$, $G_k$ is a submatrix of $G$ that contains the contribution of only sources in region $P_k$, and $F_{S_k}$ represents the log-partition function associated to the region $k$, and is defined as:

$$F_{S_k}(G_k^t \tilde{\lambda}) = \log \int \exp(\lambda^t G_k r_k) \mu(r_k | S_k)$$

$$\tilde{\pi}(S) = \frac{\prod_{k=1}^{K} \tilde{\pi}(S_k) \exp F_{S_k}(G_k^t \tilde{\lambda})}{\sum_{S} \prod_{k=1}^{K} \tilde{\pi}(S_k) \exp F_{S_k}(G_k^t \tilde{\lambda})}$$

For independent regions, the posterior conditional probability equals:

$$\tilde{\pi}(S_k) = \frac{\pi(S_k) \exp F_{S_k}(G_k^t \tilde{\lambda})}{\sum_{S_k} \pi(S_k) \exp F_{S_k}(G_k^t \tilde{\lambda})} \quad (7)$$

The joint reference law is considered as:

$$\begin{cases} \mu(r_k|S_k = 0) = \delta(r_k) \\ \mu(r_k|S_k = 1) = N(\nu_k, \Sigma_k) \end{cases} \quad (8)$$

where $\delta$ refers to the Dirac function, and $N(\nu_k, \Sigma_k)$ is a Gaussian distribution that explains the dipole's intensities inside an active patch. Let's $\alpha_k$ be the probability of the region $P_k$ to be active, i.e. $\pi(S_k = 1) = \alpha_k$. We chose zero mean ($\nu_k = 0$) and a covariance matrix ($\Sigma_k$) like in Chowdhury et al. [3] and Friston et al. [6]. The $\alpha$'s are initialized using the mean of the multivariate source pre-localization [7] scores inside the parcels obtained by the dMRI parcellation and is then updated using (7).

## 3. IMAGE ACQUISITION AND PROCESSING

Structural and diffusion MRI data were taken from 11 healthy subjects [8, 9]. The T1 weighted images of size 256 × 256 × 192 were acquired by Siemens 3T Trio with GRAPPA 3D MPRAGE sequence (TR = 2250 ms; TE = 2.99 ms; flip-angle = 9 degree; acceleration factor = 2) at 1 mm isotropic resolution. The diffusion weighted images of size 96 × 96 × 68 were collected by the same scanner at 2 mm isotropic resolution (64 gradient directions and b-value =1000 s/mm²), with one b0 image.

The W/G matter interfaces were extracted using Freesurfer [10] from T1 and remeshed to 10⁴ vertices using Brainstorm [11], then projected from the anatomical to the diffusion space. The projection was obtained by registering the brain in the two spaces using FSL [12]. The MEG/EEG forward problem, based on Boundary Element Method, is obtained using OpenMEEG [13] with a brain-skull conductivity ratio set to 80.

MEG and EEG were measured in magnetically shielded room using Elekta Neuromag Vectorview 306 system (102 magnetometers, 204 planar gradiometers), a 70 easy channel Easycap EEG cap was used to record the EEG data simultaneously. Stimuli were presented in six 7.5 min runs at 1100 Hz sampling rate. The face stimuli contains two sets of 300 gray-scale photographs, half from unfamiliar people (unknown to the subjects) and the remaining from famous people. For the third condition, 150 photographs of scrambled faces are obtained from either famous or unfamiliar people.

Across subjects, between 880 and 889 epochs were extracted (around 295 per condition i.e. Famous (Fs), Unfamiliar (Ur), Scrambled (Sd)) [8, 9]. For each of the 11 participants and conditions, we averaged the EEG and MEG measurements. This allows us to increase the signal-to-noise ratio
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** The cross-correlation matrix (CM) between all the connectivity profiles of the sources in the left ST lobe. A column or raw *i* (symmetric matrix) corresponds to the correlation between the connectivity profile of source *i* and all the others in the pre-parcel.</figcaption></figure>

and have an average overview on the cortical regions responsible of the face recognition. More details of the paradigm can be found on [9]. Low-pass filter from MNE-Python package [14] with a cutoff frequency equal to 45 Hz was used. Due to the page limitation, we show the results of only 4 subjects (2 male, 2 female, mean age 26.75).

### 3.1. Cortex parcellation

Due to the speed and memory limitations caused by the high number of sources (10⁴), we use pre-parcellated cortical surface as an initial parcellation. We use the same approach as described in Anwander et al. [4], but instead of using Brodmann atlas as a pre-parcellation, we use Destrieux atlas provided by FreeSurfer. For each Destrieux region, we compute the connectivity profiles (tractograms) of the cortical sources (seeds).

The connectivity profile of a seed *i*, a vector of size *D*, is the connectivity strength, i.e. the number of samples drawn from *i* that arrives to a voxel *j*, between seed *i* and all of the other voxels in the image space. It represents the distribution of the fibers that starts from *i*. It is obtained via a probabilistic tractography using FDT [15] toolbox from FSL. We define

<figure><img src="image_3.png" /><figcaption>**Fig. 2:** The resulting parcellation of Subject 1 with $C_{th} = 65$.</figcaption></figure>

dMRI patch as a region in which sources have similar connectivity profiles. We apply *k*-means on the cross-correlation matrix (CM) to parcellate each of the pre-parcels. The number of regions (clusters), $n_i$, inside the *i*-th pre-parcels is chosen so that the $n_i$ first biggest eigenvalues of CM represent more than a threshold $C_{th}$ (in %) of the total energy of CM.

<table><caption>**Table 1:** The total number of regions for different subjects and $C_{th}$ values. Close numbers across subjects.</caption><thead><tr><th>Subject</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th scope="row">$C_{th}$</th><td></td><td></td><td></td><td></td></tr><tr><th scope="row">60</th><td>308</td><td>297</td><td>295</td><td>304</td></tr><tr><th scope="row">70</th><td>431</td><td>400</td><td>415</td><td>432</td></tr><tr><th scope="row">80</th><td>503</td><td>544</td><td>537</td><td>533</td></tr></tbody></table>

This approach does not guarantee the regions to be spatially connected. Fig. 1 shows the correlation between pairs of connectivity profiles inside the left Superior Temporal (ST) lobe, ST can be divided into sub-regions with similar connectivity profiles. We follow the same procedure with the rest of the pre-parcels. Table 1 shows the resulting number of parcels for four different subjects and different threshold values. As expected, the number increases with the threshold value, for $C_{th} = 100$ every region contains only one source.

### 4. RESULTS AND DISCUSSION

Fig. 3 (a) compares the source reconstruction of the MEM and Minimum norm estimate (MNE) [1] on synthetic data in which one region in the Precentral gyrus (2.47 $cm^2$) was activated. The ground truth is a *sine* wave. In Fig. 3 (b), we show the effect of an additive Gaussian noise on the source reconstruction of both approaches by computing the mean squared error (MSE) between the ground truth and the estimates.

<figure><img src="image_4.png" /><figcaption>**Fig. 3:** Results of synthetic data, (a) the source reconstruction using MNE (red) and MEM (green) and (b) the MSE versus SNR for both MNE (blue) and MEM (green). MNE gives an average error of around 10% whereas MEM decreases until less than 0.5% error.</figcaption></figure>

The MEM provides more accurate results and is less affected by noise than the MN estimate. The regularization in MNE is the $l_2$ norm on the source (the trade off parameter is fixed by the Generalized Cross Validation). This smooths the intensities which reduces the reconstructed values.

The framework is tested now using the real data presented in section 3. $C_{th}$ is set to 70 for all subject. Fig. 4 shows the entropy drop, equation (1), over time for 2 subjects and different conditions. Fs-Ur represents the measurements obtained for famous photographs subtracted to the ones obtained
<figure><img src="image_5.png" /><figcaption>Fig. 4: Entropy drop due to the MEG measurements of different classes and subjects. The decrease in the entropy drop for the Fs, Ur, Sd starts around 70 ms.</figcaption></figure>

for unfamiliar photographs, the same for Fs-Sd, and Ur-Sd. The entropy curves starts significantly decreasing at around 70 ms which coincide with the literature about the visually evoked potential time. The entropy drop of Fs-Ur measurements is smaller than the other classes which means that the information brought by the MEG measurements was reduced by subtracting the measurements of Fs to Ur. Fig.5 shows the averaged absolute dipoles intensities (AADI) for Fs-Sd of the different subjects between 60 ms and 400 ms after the stimulus.

$$AADI = \frac{\sum_{t=60\text{ms}}^{400\text{ms}} abs(r_t)}{\text{number of time samples}}$$

This measure gives us the average activation of dipoles over the time window of interest and gives us an insight about which cortical regions are activated during face recognition task.

Sources in Fig.5 (c) and (d) are obtained using the Minimum Norm estimate [1] on MEG and EEG measurements respectively. The MEM-dMRI framework is used to reconstruct sources intensity from the MEG (Fig.5 (a)) and EEG measurements (Fig.5 (b)) of 4 subjects. MNE smooths the sources over the cortex which make it hard to find the real activated regions, whereas MEM provided a focal solution. The later solution is consistent with previous studies and results concerning face recognition network as shown in [16], and also to the fMRI and MEG source reconstruction of the same data set in Henson et al. [8] but with different source reconstruction algorithm in which the Inferior Occipital Gyrus (IOG) and the Fusiform gyrus (FFG) region (see Fig. 5 (e)) were found to be active. Some studies reported bilateral activation [17, 18], Subject 1-3. Others reported right sided prefrontal activation [19], see the EEG reconstruction of Subject 4 (panel 4 in Fig. 5 (b)).

Although both MEG and EEG measures the same neurophysiological processes, they are different. MEG is less affected by Skull and Scalp than EEG. The latter is sensitive to the tangential and radial sources, whereas MEG is sensitive to tangential sources. These lead to differences in the reconstruction between the MEG (Fig. 5 (a)) and EEG (Fig. 5 (b)).

<figure><img src="image_6.png" /><figcaption>Fig. 5: Averaged absolute dipoles intensities estimated by: (a) MEM on MEG, (b) MEM on EEG, (c) MNE on MEG, (d) MNE on EEG for different subjects. All values were normalized to be between 0 and 1. Numbers 1-4 represents the subject number, (e) location of IOG in orange and FFG in blue in the cortex. (f) the different views in each box.</figcaption></figure>
## 5. CONCLUSION

This paper presents a way to fuse a high spatial resolution imaging (dMRI) information to the MEM framework. This allows us to define one parcellation per subject and to use the focality advantage of the MEM. This is useful to compare the activated regions between different tasks and for multi-subject studies. Sources inside each region obtained from the dMRI informed parcellisation of the cortical substrate are assumed to follow Normal distribution and functionally independent. Synthetic and real data were used to test the reconstruction algorithm and compared with the reconstruction of MNE. For real MEG/EEG data, the results of the framework was found to be consistent with what we can find in the literature about face recognition task. For future work, MEM-dMRI will be used for functionally dependent regions case. The effect of adding a temporal regularization should be investigated.

## 6. ACKNOWLEDGMENT

We acknowledge Alexandre Gramfort for providing the code used for EEG/MEG pre-processing.

## REFERENCES

[1] R. Grech, T. Cassar, et al., “Review on solving the inverse problem in EEG source analysis,” *Journal of NeuroEngineering and Rehabilitation*, vol. 5, no. 1, pp. 25, 2008.

[2] B. Cottereau, K. Jerbi, et al., “Multiresolution imaging of MEG cortical sources using an explicit piecewise model,” *NeuroImage*, vol. 38, no. 3, pp. 439–451, 2007.

[3] R. A. Chowdhury, J. M. Lina, et al., “MEG source localization of spatially extended generators of epileptic activity: Comparing entropic and hierarchical bayesian approaches,” *PLoS ONE*, vol. 8, no. 2, pp. e55969, 2013.

[4] A. Anwander, M. Tittgemeyer, et al., “Connectivity-based parcellation of Broca’s area,” *Cerebral Cortex*, vol. 17, no. 4, pp. 816–825, 2007.

[5] C. Amblard, E. Lapalme, et al., “Biomagnetic source detection by maximum entropy and graphical models,” *Biomedical Engineering, IEEE Transactions on*, vol. 51, no. 3, pp. 427–442, 2004.

[6] K. Friston, L Harrison, et al., “Multiple sparse priors for the M/EEG inverse problem,” *NeuroImage*, vol. 39, no. 3, pp. 1104–1120, 2008.

[7] J. Mattout, M. Pélégrini-Issac, et al., “Multivariate source prelocalization (MSP): Use of functionally informed basis functions for better conditioning the MEG inverse problem,” *NeuroImage*, vol. 26, no. 2, pp. 356–373, 2005.

[8] R. N. Henson, D. G. Wakeman, et al., “A parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration,” *Frontiers in Human Neuroscience*, vol. 5, no. 76, 2011.

[9] D. G. Wakeman and R. N. Henson, “A multi-subject, multi-modal human neuroimaging dataset,” *Scientific Data*, 2015.

[10] A.M. Dale, B. Fischl, et al., “Cortical surface-based analysis I: Segmentation and surface reconstruction,” *NeuroImage*, vol. 9, pp. 179194, 1999.

[11] F. Tadel, S. Baillet, et al., “Brainstorm: A user-friendly application for MEG/EEG analysis,” *Computational Intelligence and Neuroscience*, vol. 2011, no. 1, pp. 19, 2011.

[12] S. M. Smith, M. Jenkinson, et al., “Advances in functional and structural MR image analysis and implementation as FSL,” *NeuroImage*, vol. 23, Supplement 1, no. 0, pp. S208–S219, 2004, Mathematics in Brain Imaging.

[13] A. Gramfort, T. Papadopoulo, et al., “OpenMEEG: opensource software for quasistatic bioelectromagnetics,” *BioMedical Engineering OnLine*, vol. 9, no. 1, pp. 45, 2010.

[14] A. Gramfort, M. Luessi, et al., “MEG and EEG data analysis with mne-python,” *Frontiers in Neuroscience*, vol. 7, no. 267, 2013.

[15] T.E.J. Behrens, H. Johansen Berg, et al., “Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?,” *NeuroImage*, vol. 34, no. 1, pp. 144–155, 2007.

[16] B. Rossion, R. Caldara, et al., “A network of occipito-temporal face-sensitive areas besides the right middle fusiform gyrus is necessary for normal face processing,” *Brain*, vol. 126, no. 11, pp. 2381–2395, 2003.

[17] A. Ishai and E. Yago, “Recognition memory of newly learned faces,” *Brain Research Bulletin*, vol. 71, no. 13, pp. 167–173, 2006.

[18] J.V. Haxby, L.G. Ungerleider, et al., “The effect of face inversion on activity in human neural systems for face and object perception,” *Neuron*, vol. 22, no. 1, pp. 189–199, 1999.

[19] D.J. Bayle and M.J. Taylor, “Attention inhibition of early cortical activation to fearful faces,” *Brain Research*, vol. 1313, no. 0, pp. 113–123, 2010.