Multivariate autoregressive model constrained by anatomical connectivity to reconstruct focal sources

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/mm2 ), 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 3 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) ≥ 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 ∈ R n×T ) are related to the brain sources (J ∈ R m×T ) by the following linear relationship; where n and m are the number of sensors and sources respectively, T is the length of the time window, is Gaussian observation noise and G ∈ R n×m is the lead field matrix.We assume that the sources (J t ) at time t are related to their past values by: where A is the MAR matrix and ω 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; where U is some prior on sources and λ is a positive regularization parameter.In focal activation, most of the sources are inactive.This sparsity can be obtained by putting (J) = J 21 , where J 21 = i 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 mixednorm 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: 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; where A = vec(A) and J 1, T = vec(J 1, T ) and , where 0 T ×m is a T × m zero matrix and (• ) is the transpose operator.The solution of Eq. ( 2) is; Because of the block of zeros in J * , it is possible to estimate separately each row A i of A: ., 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 AJ 0, T −1 in the data fit term of Eq. ( 1) and U (J) by U (J 0, T −1 ).Eq. ( 1) becomes: 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), 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 − GAJ 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.JT ← Ã JT −1

13:
end for 14: return J, Ã 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.
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, ±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 ±5 dB.
The reconstruction error between the estimated sources and the ground truth, J − J  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  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.

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.

Fig. 1 :
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.
(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 ≥ 5 dB (Fig. 2 (k)) and it underestimates the MAR coefficients for 5 ≥ SNR ≥ -5 dB (Fig. 2 (l)) because of the noise but they are still close to the ground truth.

Fig. 2 :
Fig. 2: The source reconstruction using mixed-norm estimates at SNR equal to (a) 40 dB, (b) 5dB and (c) -5dB.The estimated eigenvalues of Ã 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.

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.