# Multivariate Autoregressive Model Constrained by Anatomical Connectivity to Reconstruct Focal Sources

Brahim Belaoucha, Mouloud Kachouane, Théodore Papadopoulos

► To cite this version:

Brahim Belaoucha, Mouloud Kachouane, Théodore Papadopoulos. Multivariate Autoregressive Model Constrained by Anatomical Connectivity to Reconstruct Focal Sources. 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), Aug 2016, Orlando, United States. hal-01357167

HAL Id: hal-01357167

https://hal.science/hal-01357167

Submitted on 29 Aug 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
# Multivariate Autoregressive Model Constrained by Anatomical Connectivity to Reconstruct Focal Sources

Brahim Belaoucha<fnref n="1" />, Mouloud Kachouane<fnref n="1" />,<fnref n="2" /> and Théodore Papadopoulos<fnref n="1" />

**Abstract**— In this paper, we present a framework to reconstruct spatially localized sources from Magnetoencephalography (MEG)/Electroencephalography (EEG) using spatiotemporal constraint. The source dynamics are represented by a Multivariate Autoregressive (MAR) model whose matrix elements are constrained by the anatomical connectivity obtained from diffusion Magnetic Resonance Imaging (dMRI). The framework assumes that the whole brain dynamic follows a constant MAR model in a time window of interest. The source activations and the MAR model parameters are estimated iteratively. We could confirm the accuracy of the framework using simulation experiments in both high and low noise levels. The proposed framework outperforms the two-stage approach.

## I. INTRODUCTION

EEG and MEG are two non-invasive modalities that allow us to measure brain activity with high temporal resolution. Obtaining the distributed sources activation from these measurements is an under-determined problem due to the small number of measurements with respect to the number of sources. To obtain a unique solution, different constraints can be added [1]. These can be divided into three different categories: spatial, temporal, and spatiotemporal constraints.

There are several studies that show that the reconstruction improves when including spatiotemporal constraints [2], [3]. In this work, we assume that the source signals follow a multivariate autoregressive (MAR) model. Fukushima et al. [3] proposed a MAR model in which matrix elements depend on the anatomical connectivity, but their method is time-consuming and hard to implement. Other studies, [4], [5], assume that sources increase or decrease by the same factor or interact with only their direct neighbors. These assumptions significantly reduce the degrees of freedom of the dynamics. We provide more degrees of freedom by allowing source interactions constrained by anatomical connectivity obtained from dMRI information. This allows us to generalize our work to the whole brain. Focal reconstructions are important because they are consistent with the notion of functional specialization of small cortical regions. Focal Underdetermined System Solver (FOCUSS) [6] is used to detect spatially localized sources. It is applied to each time sample and tends to set to zero sources with low activation at each time sample. This makes it unsuitable for MAR parameter estimation. Our framework estimates the source activation and the MAR parameters iteratively. This approach corrects the false activation and functional connections that

can be obtained by two-stage approaches when working with low signal-to-noise ratio (SNR).

We have developed a framework constrained by a temporal dynamic over the whole brain to reconstruct sources. The source dynamics represented by the MAR model is based on the underlying anatomical connections obtained from dMRI information. Through simulations at different noise levels, we compare the accuracy of the proposed framework in the source reconstruction and MAR model parameters estimation with the ones obtained from a non-dynamical sparse solution based on mixed-norm estimates (MxNE) [7]. Our framework uses the MxNE estimates to initialize the source intensities then estimates the MAR entries iteratively.

## II. DATA ACQUISITION AND PROCESSING

Structural and diffusion MRI data were taken from Wakeman et al. [8]. Cambridge University Psychological Ethics Committee approved all experimental procedures involving human subjects [8]. The T1 weighted images of size 256 × 256×192 were acquired by a Siemens 3T Trio with GRAPPA 3D MPRAGE sequence (TR = 2250 ms; TE = 2.99 ms; flip-angle = 9°; 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<fnref n="2" />), with one b0 image. The cortical surface was extracted using Freesurfer [9] from T1 and remeshed to 2003 vertices (sources). The transformation between the anatomical and diffusion space was obtained by registering the brain in the two spaces using FSL [10]. The MEG/EEG forward problem, based on Boundary Element Method (BEM), is obtained using OpenMEEG [11], [12].

## III. STRUCTURAL CONNECTIVITY

The structural connectivity of the sources (nodes on the mesh) is obtained from the dMRI information. Each node in the cortical mesh is considered as a seed. For each seed, we compute the connectivity profile, i.e. its connectivity with every voxel in the diffusion image space. Connectivity profiles are obtained by running probabilistic tractography for all seeds. We assume a model of multiple fiber orientations for each voxel developed in FSL [14] with 5 × 10³ samples. The connectivity between node *i* and *j*, $c_{ij}$, is the number of samples that were drawn from *i* and arrived at *j* divided by the total number of samples. *CM* is the connectivity matrix with elements $CM(i,j) = \max(c_{ij}, c_{ji})$. The source *i* and *j* are considered to be structurally connected if $CM(i,j) \ge$

<footnote n="1">¹ Université Côte d'Azur, Inria, Nice, France belaoucha.brahim@etu.unice.fr brahim.belaoucha@inria.fr</footnote>

<footnote n="2">² USTHB, LRPE, ParIMed, Algiers, Algeria</footnote>
$c_{th}$. We set $c_{th}$ to 0.1 to remove weak connections. We obtained 3935 anatomical connections.

IV. ALGORITHM

The EEG/MEG measurements ($M \in \mathbb{R}^{n \times T}$) are related to the brain sources ($J \in \mathbb{R}^{m \times T}$) by the following linear relationship;

$$M = GJ + \epsilon,$$

where $n$ and $m$ are the number of sensors and sources respectively, $T$ is the length of the time window, $\epsilon$ is Gaussian observation noise and $G \in \mathbb{R}^{n \times m}$ is the lead field matrix. We assume that the sources ($J_t$) at time $t$ are related to their past values by:

$$J_t = AJ_{t-1} + \omega_t,$$

where $A$ is the MAR matrix and $\omega_t$ is Gaussian generator noise. The elements of the MAR matrix are constrained by the structural connectivity obtained from dMRI information. The only non-zero elements of $A$ represent anatomically connected source pairs including self-connections. The dynamics of the model depend on the eigenvalues of $A$. The sources can be obtained by minimizing the following functional;

$$E(J) = \frac{1}{2T} \|M - GJ\|_2^2 + \lambda U(J), \quad (1)$$

where $U$ is some prior on sources and $\lambda$ is a positive regularization parameter. In focal activation, most of the sources are inactive. This sparsity can be obtained by putting $U(J) = \|J\|_{21}$, where $\|J\|_{21} = \sum_i \sqrt{\sum_t J_t(i)^2}$. This norm, $l_{12}$, induces sparsity over the time window i.e. the sites of active sources will be the same over the entire time window of interest. The solution of Eq. (1) is called the mixed-norm estimate (MxNE). More details about the MxNE can be found in [7]. To enhance the solution of the mixed-norm solution in low SNR values and to estimate the functional connectivity which is represented by the MAR model, we introduce an iterative framework. The MAR matrix, $A$, is estimated over a time window of interest by minimizing $D$:

$$D(A) = \|J_{1,T} - AJ_{0,T-1}\|_2^2,$$

where $J_{1,T}$ and $J_{0,T-1}$ are the values of the sources in a time window of $T$ samples. This can be rewritten as;

$$D(\mathbf{A}) = \|\mathbf{J}_{1,T} - \mathbf{J}^*\mathbf{A}\|_2^2, \quad (2)$$

where $\mathbf{A} = \text{vec}(A)$ and $\mathbf{J}_{1,T} = \text{vec}(J_{1,T})$ and

$$J^* = \begin{pmatrix} J'_{0,T-1} & 0_{T \times m} & \cdots & 0_{T \times m} \\ 0_{T \times m} & J'_{0,T-1} & \cdots & 0_{T \times m} \\ \vdots & \vdots & \ddots & \vdots \\ 0_{T \times m} & 0_{T \times m} & \cdots & J'_{0,T-1} \end{pmatrix},$$

where $0_{T \times m}$ is a $T \times m$ zero matrix and ($'$) is the transpose operator. The solution of Eq. (2) is;

$$\mathbf{A} = (J^{*'} J^{*})^{-1} J^{*'} \mathbf{J}_{1,T}, \quad (3)$$

Because of the block of zeros in $J^*$, it is possible to estimate separately each row $A_i$ of $A$:

$A'_i = (J_{0,T-1} J'_{0,T-1})^{-1} J_{0,T-1} J_{1,T}(i) \quad i \in \{1, \dots, m\}$, where $J_{1,T}(i)$ is the $i^{th}$ row of $J_{1,T}$. To estimate the source intensities and their interaction iteratively, we replace $J_{1,T}$ by $A J_{0,T-1}$ in the data fit term of Eq. (1) and $U(J)$ by $U(J_{0,T-1})$. Eq. (1) becomes:

$$E(J_{0,T-1}) = \frac{1}{2T} \|M_{1,T} - GA J_{0,T-1}\|_2^2 + \lambda \|J_{0,T-1}\|_{21} \quad (4)$$

The framework, iterative Source and Dynamics Reconstruction (iSDR), is explained in Algorithm 1. We start by initializing the MAR matrix $A$ with the identity $(I)$, $\tilde{J}$ with the MxNE solution and the maximum number of iterations $N$. In step 4, we estimate the sources by optimizing Eq. (1). Then $A$ is updated by Eq. (3) in steps 6-8. The algorithm stops iterating if there is a small change in the data fit error, $RSS = \|M_{1,T} - GA J_{0,T-1}\|_2$, or after $N$ iterations. The difference between Eq. (4) and Eq. (1) (MxNE) is that the lead field matrix, $G$, is multiplied by the MAR matrix, $A$. The columns and rows of $A$ corresponding to inactive sources are set to zero in step 7. This allows the framework to reduce the effect of inactive sources on the measurements, i.e. reduce the number of unknowns resulting in better estimation of the source activation.

**Algorithm 1** iterative Source and Dynamics Reconstruction (iSDR)

1: **procedure** iSDR($G, M, N$)

2: **Initialize:**
$\tilde{A} \leftarrow I$, $(n, m) \leftarrow \text{shape}(G)$, $\varepsilon$

3: **for** $i=1:N$ **do**

4: $\tilde{J}_{0,T-1} \leftarrow \arg \min_{E} (J_{0,T-1}(E)) \triangleright \text{Estimate } \tilde{J}$

5: $\tilde{J}_T \leftarrow \tilde{A} \tilde{J}_{T-1}$

6: **for** $j=1:m$ **do**
$\tilde{A}'_j \leftarrow (J_{0,T-1} J'_{0,T-1})^{-1} J_{0,T-1} J_{1,T}(j)$
$\triangleright \text{Update } \tilde{A}$

7: $\tilde{A}'_j \leftarrow (J_{0,T-1} J'_{0,T-1})^{-1} J_{0,T-1} J_{1,T}(j)$

8: **end for**

9: $RSS^i \leftarrow \|M_{1,T} - G\tilde{A}\tilde{J}_i\|_2$
10: **if** $RSS^i - RSS^{i-1} < \varepsilon$ **then**

11: **break**

12: **end if**

13: **end for**

14: **return** $\tilde{J}$, $\tilde{A}$

15: **end procedure**

V. RESULTS AND DISCUSSION

We used a realistic three layers head model. 70 EEG electrodes were considered. To test the framework, we first obtain the structural connectivity from the dMRI as explained in section III. Then, we investigate the parameter estimation accuracy of the proposed framework by using simulated data generated from a simple MAR model. Four sources ($S1$, $S2$, $S3$, and $S4$) are located in the Occipital lobe, see Fig. 1 (a) and connected anatomically, see Fig. 1 (b). We activated $S1$ and $S2$ while $S3$, $S4$, and the other sources are inactive. We decided not to activate $S3$ and $S4$ to investigate the effect of the assumption that the sources functionally interact if they are anatomically connected. $S1$ and $S3$ are anatomically connected but functionally independent.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Location of the simulated sources on the inflated cortical surface (a) and their anatomical connections (b), only the connections between the four sources. The magnitude of the simulated sources (c), sampling frequency equal to 1 KHz.</figcaption></figure>

$$ \begin{cases} J_t(S1) = 0.96J_{t-1}(S1) + 0.25J_{t-1}(S2) \\ J_t(S2) = -0.25J_{t-1}(S1) + 0.95J_{t-1}(S2) \end{cases} $$

with $J_0(S1) = -J_0(S2) = 5$ nAm (nano Ampere), the simulated sources can be seen in Fig. 1 (c). The MAR matrix has two non-zero eigenvalues $0.955 \pm 0.25j$. The regularization parameter, $\lambda$, is fixed by a 3-fold cross validation [13].

We simulate sources in low and high SNR. The simulated measurements were obtained from the simulated sources (Fig. 1 (c)) multiplied by the lead field matrix, $G$, and we added a Gaussian observation noise [SNR: 40, $\pm$5 dB] of zero mean. The source reconstruction of MxNE can be seen in Fig. 2 (a) to (c). Fig. 2 (d) shows the estimated eigenvalues using two-stage approach (MxNE then Eq. (3)) at SNR equal to 5 dB. Fig. 2 (e) to (g) represent the result of the proposed method at different noise levels.

Fig. 2 (h) to (j) show the true, in black, and the estimated, in red, eigenvalues by the proposed framework of the MAR matrix in polar coordinates at different noise levels and for 100 runs. The majority of the eigenvalues are estimated to be zero, this can be guessed clearly from the reconstructions of the iterative approach at the different noise levels, Fig. 2 (e) to (g), in which the framework detected the true inactive sources (the corresponding MAR coefficients equal to zero). Fig. 2 (k) and (l) represent the mean over 100 runs of the estimated functional interactions between the four sources at SNR equal to $\pm$5 dB.

The reconstruction error between the estimated sources and the ground truth, $\|\tilde{J} - J\|_2^2$, of the MxNE estimates and the proposed framework (iSDR) at different noise levels is

<table><caption>TABLE I: The mean and standard deviation (std) of the reconstruction error over 100 runs at different noise levels (SNR = {40,5,-5} dB). The values are in nAm.</caption><thead><tr><th rowspan="2"></th><th colspan="2">40 dB</th><th colspan="2">5 dB</th><th colspan="2">-5 dB</th></tr><tr><th>mean</th><th>std</th><th>mean</th><th>std</th><th>mean</th><th>std</th></tr></thead><tbody><tr><td>MxNE</td><td>11.19</td><td>0.008</td><td>38.45</td><td>0.47</td><td>237</td><td>2.59</td></tr><tr><td>iSDR</td><td>0.15</td><td>0.0045</td><td>5.01</td><td>0.19</td><td>15.31</td><td>0.61</td></tr></tbody></table>

shown in Table 1. The mean error is computed over 100 runs at the different noise levels. The iSDR enhances significantly the mixed-norm estimates solution especially at low SNR settings, compare Fig. 2 (b) with (f) and (c) with (g) and see the last column of Table 1.

The iSDR provides accurate results because MxNE sets sources to zero in the whole time window of interest. This results in rows and columns of zeros in the MAR matrix, step 7 in Algorithm. In the next iteration, $A$, with the zero columns and rows, is multiplied by $G$. This is equivalent to removing some generators from the source space resulting in a better approximation. The number of iterations, $N$, is inversely proportional to noise level i.e. high $N$ value for low SNR.

The two-stage approach found false-positive active sources in low SNR, Fig. 2 (b) and (c), which result in false functional connections. This is the limitation of two-stage approximation of the MAR model. The iterative framework could correct the false activation, see Fig. 2 (f) and (g). It could detect the correct connections between the sources for SNR $\ge$ 5 dB (Fig. 2 (k)) and it underestimates the MAR coefficients for $5 \ge$ SNR $\ge$ -5 dB (Fig. 2 (l)) because of the noise but they are still close to the ground truth.

In this paper, we present an iterative framework for EEG and MEG source reconstruction constrained by a whole brain dynamics model based on dMRI information. In our simulation study, the proposed framework correctly reconstructed sources and estimated the MAR elements, while the non-dynamical two-stage approach based on MxNE solution detects false-positive sources resulting in false functional connections. The framework assumes that the source dynamics is constant over a time window of $T$ samples. The optimal number of samples needed to estimate the MAR entries is proportional to the maximal number of anatomical connections per source. It is hard for the framework to detect sources that have small activation in time and duration. Like all initialization dependent approaches, the initialization has to be within some proximity to the true sources, in order to find the correct activation. This means that as long as the initial active set contains true active sources, we are able to recover correct activations and interactions.

Representing the source dynamics by linear autoregressive modeling would introduce errors to the estimated functional interactions for sources with non-linear relationships. In this work, we dealt with focal sources. Extended generators (constant activation per region) can be assumed on regions obtained from dMRI-based parcellation [3], [14] in the proposed framework. Higher-order MAR model can be assumed.
<figure><img src="image_3.png" /><figcaption>Fig. 2: The source reconstruction using mixed-norm estimates at SNR equal to (a) 40 dB, (b) 5 dB and (c) -5 dB. The estimated eigenvalues of $\hat{A}$ using two-stage approach, i.e MxNE then Eq. (3) in (d). The source reconstruction obtained by the proposed method (iSDR) at SNR (e) 40 dB result after 4 iterations, (f) 5 dB result after 10 iterations and (g) -5 dB result after 14 iterations. The estimated eigenvalues (100 runs) by the proposed iterative method at different noise levels, (h) 40 dB, (i) 5 dB and (j) -5 dB. In (k) and (l), we show the mean over 100 runs of the estimated functional interactions between $S1$, $S2$, $S3$ and $S4$ for SNR equal to 5 and -5 dB respectively.</figcaption></figure>

## VI. CONCLUSION

We developed a source reconstruction framework constrained by a whole brain dynamical model based on dMRI information. The framework iteratively solves the source activation and MAR entries. In our study, the framework could correct the false activation in low SNR level generated by the two-stage approach and detect the correct functional interactions between sources. Application to real datasets and comparison of the solutions to the physiological literature is essential to validate this framework.

## REFERENCES

[1] R. Grech et al., "Review on solving the inverse problem in EEG source analysis," *Journal of NeuroEngineering and Rehabilitation*, vol. 5, no. 1, pp. 1-33, 2008.

[2] C. Lamus et al., "A spatiotemporal dynamic distributed solution to the MEG inverse problem," *NeuroImage*, vol. 63, no. 2, pp. 894-909, 2012.

[3] M. Fukushima et al., "MEG source reconstruction constrained by diffusion MRI based whole brain dynamical model," in *2013 IEEE 10th International Symposium on Biomedical Imaging*, April 2013, pp. 1002-1005.

[4] E. Giraldo et al., "Estimation of dynamic neural activity using a kalman filter approach based on physiological models," in *2010 Annual International Conference of the IEEE Engineering in Medicine and Biology*, Aug 2010, pp. 2914-2917.

[5] C. J. Long et al., "Large scale kalman filtering solutions to the electrophysiological source localization problem- a meg case study," in *2006 Annual International Conference of the IEEE Engineering in Medicine and Biology*, Aug 2006, pp. 4532-4535.

[6] I. F. Gorodnitsky et al., "Neuromagnetic source imaging with FO-CUSS: a recursive weighted minimum norm algorithm," *Electroencephalography and Clinical Neurophysiology*, vol. 95, no. 4, pp. 231-251, 1995.

[7] A. Gramfort et al., "Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods," *Physics in Medicine and Biology*, vol. 57, no. 7, pp. 1937-1961, 2012.

[8] D. G. Wakeman et al., "A multi-subject, multi-modal human neuroimaging dataset," *Scientific Data*, 2015.

[9] A. M. Dale et al., "Cortical surface-based analysis I: Segmentation and surface reconstruction," *NeuroImage*, vol. 9, pp. 179-194, 1999.

[10] S. M. Smith et al., "Advances in functional and structural MR image analysis and implementation as FSL," *NeuroImage*, vol. 23, Supplement 1, pp. S208-S219, 2004, Mathematics in Brain Imaging.

[11] A. Gramfort et al., "OpenMEEG: opensource software for quasistatic bioelectromagnetics," *BioMedical Engineering OnLine*, vol. 9, no. 1, pp. 1-20, 2010.

[12] J. Kybic et al., "A common formalism for the integral formulations of the forward EEG problem," *Medical Imaging, IEEE Transactions on*, vol. 24, no. 1, pp. 12-28, Jan 2005.

[13] F. Pedregosa et al., "Scikit-learn: Machine learning in Python," *Journal of Machine Learning Research*, vol. 12, pp. 2825-2830, 2011.

[14] B. Belaoucha et al., "Cortical surface parcellation via dMRI using mutual nearest neighbor condition," in *2016 IEEE International Symposium on Biomedical Imaging*, Prague, Czech Republic, Apr. 2016, pp. 903-906.