Large brain effective network from EEG/MEG data and dMR information

Over the past 30 years, neuroimaging has become a predominant technique. One might envision that over the next years it will play a major role in disclosing the brain’s functional interactions. In this work, we use information coming from diffusion magnetic resonance imaging (dMRI) to reconstruct effective brain network from two functional modalities: electroencephalography (EEG) and magnetoencephalography (MEG).


I. INTRODUCTION
MEG and EEG imaging techniques record the brain cortical activity in real time.Localizing active regions from these measurements, in the source distributed model, involves solving an ill-posed inverse problem.Prior on the source space must be set to obtain a unique solution.
Minimum norm estimate (MNE) is one of the most used inverse solvers because it leads to a linear inverse operator [1,2].It results to widely extended active regions which are, biologically, not plausible.Among other solvers, we can find the minimum current estimate (MCE).It provides a focal reconstruction but is non-smooth in time because the sparse prior is applied at each time sample independently.Mixed-norm estimate (MxNE) [3] was introduced to overcome this limitation.It provides spatial focality, by applying l 1 -norm on the source space, and temporal smoothness, by using l 2 -norm on time courses.
Functional connectivity analysis has been traditionally implemented at the sensor level, but lately, a number of studies have begun to use source-space analysis.It has the potential of providing more accurate information regarding regions' functional in-teractions [4,5,6].Some of the functional measures (e.g.Directed Transfer Function (DTF) and Partial Directed Coherence (PDC)) can be derived from multivariate autoregressive model (MAR) which is more than a combination of purely temporal and purely spatial correlation, as it also takes into account the time-lagged correlation between sources.
In [4], the authors constrain the dynamics of the sources by a MAR model.They assume that there exist only one time-lag, which depends on the length of anatomical connections, between a pair of regions.In [6], the authors use a variant of MxNE constrained with a MAR model of order 1 to reconstruct regions' interactions.In both works, complex time courses can not be reconstructed because of limits in the degree of freedom of the model.In this work, we generalize the work of [6] by constraining the dynamics of the cortical regions with higher MAR model whose matrix elements are constrained by the anatomical connectivity and neighbors.We call it iterative source and dynamics reconstruction of order p (iSDR(p)).We test and compared our method with real and synthetic data.

II. METHOD
Under the quasi-static approximation of the Maxwell's equations, the EEG/MEG data M t is a linear superposition of the source signals J t at the same time t: where G ∈ R m×n (m: size of sensor array, n: size of source space) is the lead-field matrix and ε t is N (0, Σ c ).To reduce the problem complexity, we assume a constant activation per cortical region.This reduces the size of the lead-field to m × K, where K is the number of cortical regions.We use the mutual nearest neighbor parcellation algorithm presented in [7,8] to divide the cortical surface into functional regions from dMRI.The source dynamics is constant during a time window T and it follows MAR model of order p: where ω t is N (0, Σ s ), p is a positive integer that defines the order of the MAR model, A i is a K ×K matrix which relates sources' current values to its pasts.The only nonzero effective interactions in A corresponds to either structurally connected or neighboring regions.The behavior of the MAR model is obtained from the companion matrix Φ: where I is the K × K identity matrix.The MAR model is stable (stationary) if all eigenvalues of Φ are inside the unit circle and nonstationary if it has eigenvalues on the unit circle.The model is unstable if at least one eigenvalue is outside the unit circle.
The sources estimates J are obtained by minimizing the following functional given a MAR model: where i.e. l 2 -norm over time and l 1 -norm over space.
where G i = GA i .The MAR elements are obtained by fitting the source estimates to the MAR model: We solve iteratively Eq 4 and 6 until J k − J k−1 < where J k is the reconstructed sources at iteration k.

III. DATA
MEG and EEG were measured using 306 system (102 magnetometers, 204 planar gradiometers), and 70 EEG channels recorded EEG data simultaneously.Stimuli were presented in six 7.5 min runs.The face stimuli contain two sets of 300 gray scale photographs, half from unfamiliar people (unknown to the participants) and the remaining from famous people.In a third condition, 150 photographs of scrambled faces are obtained from either famous or unfamiliar people.The reader is referred to [9] for more details.
In this work, we are interested in localizing face recognition areas.For this reason, we use only the measurement acquired when using photos of famous people subtracted to the ones obtained when using scrambled faces.Low pass filter of cut-off frequency 35 Hz was used to smooth the data.The MEG/EEG forward problem, lead field matrix G, is obtained using OpenMEEG [10,11].

A. Simulation
We use a random lead-field matrix of size 20 × 500 and whose elements are obtained from a normal distribution N (0, 1).Each region is connected randomly to 4 different regions.We activate 50 random pairs of regions, one pair at a time, 100 times at different noise levels (signal-to-noise ratio (SNR) = {15, 10, 5} dB) using the following simple MAR model: , with J 0 (S 1 ) = 6.15,J 0 (S 2 ) = −3.64, We estimate the sources and their interactions by initializing A 2 and A 1 with values to have an insight about the effect of the MAR's initialization on the resulting estimate.We use the following combinations: A 2 = βI and A 1 = αI with (α, β) ∈ {(1, 0), ( 34 , 1 4 ), ( 1 2 , 1 2 )}.The results are presented in Fig (1).Fig ( 1a) represents the reconstruction error computed as J gt − J r 2 , where J gt and J r are the simulated ground truth and reconstructed activation, using MxNE and iSDR(p=2) respectively.The estimated initial values of S 1 and S 2 obtained from iSDR(p=2) are shown in Fig ( 1b).Fig ( 1c) represents the nonzero eigenvalues of the estimated companion matrix Φ.
iSDR(p) outperforms MxNE especially in low SNR.For (α, β) = (1, 0), iSDR favors the initial values (J t=0 ) to be close to 0, see blue dots in Fig ( 1b), due to initializing A 2 = 0 i.e no interactions between t and t − 2. Setting β to nonzero improves both the temporal reconstruction of signals and their dynamics (eigenvalues of Φ). iSDR could reconstruct accurately two eigenvalues, right side of the plane in Fig ( 1c).The remaining two are affected by the choice of (α, β).The results are less dispersed when (α, β) = ( 12 , 1 2 ) (see left side of Fig ( 1c)).B. Real data iSDR(p=2) is applied to MEG data of 9 subjects from [9].In [12], the authors show that subjects may have right or bilateral fusiform gyrus (FF) activation if they are right or left handed, respectively.In Several works [9,4] show that activation in the FF between 170 and 200 ms is observed when performing face recognition task which coincides with what we observe in this work.In the majority of subjects, regions in the superior posterior temporal lobe, medial parietal, temporal pole (TP), orbitofrontal (OF) are found to be active.This coincides with the findings of other works with the same dataset [9,4].We have found interactions between TP and frontal lobe (FL), left and right FL, inferior occipital and temporal gyrus with OF.
We applied our method only with MAR model of order 2. Our approach provides functionally plausible results when tested with MEG data provided by [9].

V. CONCLUSION
In this work, we proposed a source reconstruction method constrained by multivariate autoregressive model from functional imaging EEG/MEG data.Application to real data of other mental tasks must be investigated and compared to physiological literature to validate our method.Information criterion can be used to obtain the order of MAR model.

Fig. 1 :
Fig. 1: In (a), mean, over 100 simulations, of reconstruction error Er of MxNE and iSDR(p=2) with different initialization of A and SNR levels.(b) Estimated initial values of S1 and S2 using iSDR(p = 2).(c) Eigenvalues of the companion matrix.The results are shown at 3 different SNR.Black dots represent ground truth values.

Fig 2 ,
we show the estimated temporal course of the right/left FF (R/LFF).In our work, seven subjects show higher activation in the RFF, Fig (2) left, while only two show higher activation in LFF, Fig (2) right.Due to page limitation, we show only the effective connectivity of Subject 1, Fig (3).

Fig. 3 :
Fig. 3: Effective network of A1, left, and A2, right, obtained from iSDR(p=2) and MEG data [9] during face recognition task.Red spheres represent center of active cortical regions.Edges represent uni/bidirectional interactions.Self-interactions are not shown for visualization clarity.Edge's color represents its strength.