Multi-modal EEG and fMRI Source Estimation Using Sparse Constraints

Abstract

In this paper a multi-modal approach is black presented and validated on real data to estimate the brain neuronal sources based on EEG and fMRI. Combining these two modalities can lead to source estimations with high spatio-temporal resolution. The joint method is based on the idea of linear model already presented in the literature where each of the data modalities are first modeled linearly based on the sources. Afterwards, they are integrated in a joint framework which also considers the sparsity of sources. The sources are then estimated with the proximal algorithm. The results are validated on real data and show the efficiency of the joint model compared to the uni-modal ones. We also provide a calibration solution for the system and demonstrate the effect of the parameter values for uni- and multi-modal estimations on 8 subjects.

1 Introduction

The integration of different modalities has brought up interest to understand the complex procedure of the brain. Two of these modalities are functional Magnetic Resonance Imaging (fMRI) and electroencephalography (EEG). fMRI measures the oxygenation of the blood flow, which is closely correlated to the neuronal activity. The oxygenated blood flows to the active neurons, and this leads to contrast images called Blood Oxygenation Level Dependent (BOLD). Although fMRI has a high spatial resolution, the images taken from the volume of the brain have limited temporal precision. On the other hand, EEG measures the neuronal electrical activity of the brain in the scale of milliseconds. However, it does not directly provide the spatial location of active neurons, and the localization of the active electric sources from EEG is only possible by solving an inverse problem which is ill-posed and undetermined [1]. A review of the techniques for solving the EEG inverse problem can be found in [2].

The use of both EEG and fMRI in one multi-modal framework can result in an estimation of neuronal activities with a high resolution in time and also in space. The complementary strengths that these two modalities can provide, has stimulated a variety of studies [3, 4]. These methods could be fMRI-constrained, or EEG-informed where models driven on one modality act as the spatial (or temporal) prior on the solutions of the other modality [5, 6]. EEG and fMRI can also be used in symmetrical approaches, so that the information fusion is based on both data at the same time. For example, joint Independent Component Analysis (ICA) is used to analyze both data in a joint space as in [7], or Bayesian approaches are used in this regard [8]. According to the literature, the relation of neuronal activity, hemodynamics, and fMRI is still unclear [9]; however, the existing methods model this relationship based on the balloon model or neural mass models [10]. Although the integration of these two modalities has been studied a lot, it is mainly validated according to simulated data. A symmetrical approach is proposed in [11] according to the balloon model, and the current paper adapts the proposed linear approaches to validate it on real data.

In the current study we propose to adapt a symmetrical approach based on EEG and fMRI for the estimation of the brain sources and validate it on real data. The joint algorithm is composed of two main parts: the simple linear model of both data, and the sparsity of the sources. The linear model is a first approximation and in this study we intend to show that despite the possible limits, the sparse and joint model can provide better spatio-temporal results in source reconstruction for real data, so that it can justify developing more complex models based on multi-modality.

In Sect. 2 EEG and BOLD inverse problems are explained, based on which the joint model is introduced with the model parameters. The sources are then estimated with the proximal algorithm. Section 3 explains the experiment, the multi-modal data, and its analysis, and the results based on the real data, shown in Sect. 4, validates the efficiency of the joint approach, and justifies the calibration of the model parameters. Finally the conclusions are discussed in Sect. 5.

2 Method

2.1 Model

EEG Inverse Problem. The aim is to estimate the neural sources at time instants of EEG. According to the forward problem (see Sect. 3), the $J\times L$ leadfield matrix, G, is calculated, where J and L are the number of electrodes and sources, respectively. This matrix contains coefficients that linearly link the source amplitudes to the electrode potentials. The estimation of the source matrix, $\mathbf {S}$, from the following equation, is the solution of the EEG inverse problem:

$$\begin{aligned} \mathbf {E}=\mathbf {G}\mathbf {S}+\mathbf {N}_e, \end{aligned}$$

(1)

where $\mathbf {E}$ is a $J\times T_e$ matrix corresponding to the signal obtained from J electrodes during $T_e$ instants. As noted, there is no unique solution of $\mathbf {S}$ because EEG has much fewer channels than the number of sources ($L\gg J$). That is why the hypothesis of the sparsity is made about the sources.

Fig. 1.

On the left two example of neuronal stimuli. On the right, their corresponding HR. The BOLD response is the convolution of the sources and the HRF.

BOLD Inverse Problem. In the framework of balloon model [12], the BOLD response is written as the linear convolution of the neural activity and the hemodynamic response, h(t). This is depicted in Fig. 1. Here, the electrical source locations are corresponding to voxels. So, considering $\mathbf {B}$ a $L\times T_b$ matrix as the observed BOLD signal of $T_b$ instances for L voxels, the following equation holds:

$$\begin{aligned} \mathbf {B}=\mathbf {S}\mathbf {H}+\mathbf {N}_b. \end{aligned}$$

(2)

where $\mathbf {N}_b$ is the Gaussian noise, and $\mathbf {H}$ is a $T_e\times T_b$ matrix whose $(p,q)^{th}$ element is $h(t_p-u_q)$, where $t_p$ and $u_q$ correspond to the time samples of EEG and BOLD, respectively. Note that $T_e\gg T_b$ makes the BOLD inverse problem undetermined in time, that is why the sparsity of sources is important in finding the solution.

Joint Inverse Problem. Two minimization problems can be induced from the inverse problems of Eqs. (1) and (2) to estimate sources $\mathbf {S}$, and these two can be combined symmetrically in a single model. Assuming a Gaussian noise and sparse sources we get a least square problem with $L_1$ regularization. The sparsity is a logical hypothesis since a special mental task would activate only a certain number of sources in functional brain regions. Considering the penalty term $\phi (\mathbf {S})=||\mathbf {S}||_1= \sum _{i} \sum _k\mathbf |{S}_{i,k}|$, the sources can be finally estimated as:

$$\begin{aligned} \mathbf {\hat{S}}&= \underset{\mathbf {S}}{{\text {argmin}}}\Big \{\alpha || \mathbf {G}\mathbf {S}-\mathbf {E}||_2^2+(1-\alpha )||\mathbf {S}\mathbf {H}-\mathbf {B}||_2^2+\lambda _j||\mathbf {S}||_1\Big \}, \end{aligned}$$

(3)

where $||.||_2$ denotes the $L_2$ norm, $0<\lambda _j$, and $0\le \alpha \le 1$. $\lambda _j$ adapts the solution as a good tradeoff between the joint data-fidelity term and the penalty term. $\alpha$ balances the relative importance of each modality in the joint approach.

2.2 Proximal Algorithm

Proximal algorithms are for obtaining optimization solutions for composite functions, which might not be differentiable. The joint problem has already been formulated as

$$\begin{aligned} \mathbf {\hat{S}}= \underset{\mathbf {S}}{{\text {argmin}}}\big \{f(\mathbf {S})+\lambda _j\phi ({\mathbf {S}})\big \}, \end{aligned}$$

(4)

where f is the data-fidelity term, a proper convex function whose gradient is Lipschitz continuous, while the defined $\phi$ is not differentiable. First consider the general definition of the proximity operator $\text {prox}_{\lambda ,\phi }(\mathbf {y})= \underset{\mathbf {y}}{{\text {argmin}}}\{\frac{1}{2}||\mathbf {x}-\mathbf {y}||_2^2 +\lambda \phi (\mathbf {y})\}.$ Now, an iterative method, called Iterative Soft Thresholding (IST) [13], is used with a small gradient descent step $\mu$ to find the solution as follows:

$$\begin{aligned} \mathbf {S}^{(k+1)}=\text {prox}_{\mu ,\phi }\Big (\mathbf {S}^{(k)}-\mu \nabla f(\mathbf {S}^{(k)})\Big ). \end{aligned}$$

(5)

It is proven that this algorithm converges to a minimum [14].

Other norms (i.e. $L_2$ or mixed $L_{12}$) could be used for the $\phi$ regularization to avoid the estimation of scattered sources [11]; but, comparing different penalty terms is beyond the scope of this paper. Interested readers can refer to [11] for a theoretical investigation of the differences between those terms in that context. Having the $L_1$ norm penalty, the proximity operator on the elements of the source matrix, $\mathbf {S} \in \mathbb {R}^{V\times T}$ is given by the following equation:

$$\begin{aligned}{}[\text {prox}_{\lambda ,||.||_1}(\mathbf {S})]_{v,t}=\frac{{\mathbf {S}}_{v,t}}{|\mathbf {S}_{v,t}|} (|\mathbf {S}_{v,t}|-\lambda ) \end{aligned}$$

(6)

2.3 Parameter Calibration

Consider again the joint model defined in Eq. (3). Indeed, $\alpha$ is hard to be tuned between the two sides of the data-fidelity term, since these two data-fidelity terms have no reason to have the same order of magnitude; so, a calibration of the respective roles of the two signals should be embedded. Therefore, $\lambda _e$ and $\lambda _b$ are separately defined for each data-fidelity term (corresponding to EEG and BOLD, respectively). The joint model can thus be re-written considering the new calibration parameters:

$$\begin{aligned} \mathbf {\hat{S}}= \underset{\mathbf {S}}{{\text {argmin}}}\Big \{ \frac{\alpha }{\lambda _e \lambda _j}||\mathbf {G}\mathbf {S}-\mathbf {E}||_2^2 + \frac{(1-\alpha )}{\lambda _b \lambda _j} ||\mathbf {S}\mathbf {H}-\mathbf {B}||^2_2 + \phi (\mathbf {S})\Big \} , \;\; \alpha \in [0, 1]. \end{aligned}$$

(7)

Note that if $\alpha$ is 0 (or 1), the model will be uni-modal based on BOLD (or EEG). So, if $\alpha \in (0,1)$ (ideally $\alpha =0.5$) gives better estimations than uni-modal approaches, this can validate the joint method.

3 Experimental Data and Implementations

The simultaneous EEG and fMRI are recorded from 8 healthy right-handed subjects while doing a motor task of clenching the right hand repeatedly. A study session consists of 10 blocks of 20 s from which 5 are task each followed by a rest block. EEG is obtained via 64 electrodes based on the standard 10–20 system with $fs=200$  Hz. The ECG (Electrocardiogram) and the gradient artifacts caused by MR system are separately eliminated. The sources are estimated for the first 15 s of a task block, since at the end of each block the subjects might not guarantee to do the task well because of fatigue. At the same time, the subject would enter a 3T MRI system, for the acquisition of the structural and functional MRI. Every 2 s a volume of the brain is being recorded. The fMRI data is then spatially corrected for the head motion.

The head model, is made up of 3 nested concentric spheres with homogeneous and isotropic conductivity representing the brain, the skull, and the scalp [15]. The electrodes are placed on the scalp sphere surface defined in the MRI anatomical space (according to the 10–20 position). The sources are modeled by fixed dipoles located at the voxels’ locations that are selected by brain extraction and segmentation (to consider the gray matters as the possible source locations). These sources are then down-sampled so that their number correspond to the physiological number of the current sources. A sphere, representing the brain, is then fitted to the sources. Finally, after defining another sphere for the skull, the subject-specific 3-layer spherical model is defined with different conductivities considered for each of the layers according to [15, 16]. The forward problem is then solved to linearly correlate the sources and the observed EEG potentials in the lead-field matrix.

4 Results

The subjects were involved in the task of right hand movement, so we expect the activation on the motor cortex of the left hemisphere. The motor cortex is composed of the primary motor cortex (M1), premotor cortex (called PMA), the supplementary motor area (SMA), and the posterior parietal cortex (PPC). The brain space is aligned with the AAL (Automated Anatomical Labeling) atlas, based on which the brain region in which each source is placed is determined. When the sources $\mathbf {S}$ are estimated, the power p(v) of each source v is calculated as the $L_1$ norm of its activity during all task duration: $p(v)=||\mathbf {S}(v,:)||$. So, the power score of each region r is defined as

$$\begin{aligned} \text {power}(r) = \sum _{v\in r}{p(v)}. \end{aligned}$$

(8)

Parameter Estimation. The performance, noted Q, of any estimated source matrix is defined to be higher for an estimation with a power score higher in the left motor cortex and less in other regions:

$$\begin{aligned} Q(\mathbf {S})=\sum _{r=\{M1,PMA,PPC\}_{left}} {\text {Power}(r)} -\sum _{r\ne \{M1,PMA,PPC\}_{left}} {\text {Power}(r)} \end{aligned}$$

(9)

So, the estimation and calibration of parameters is done based on the performance. $\lambda _b$ and $\lambda _e$ are first estimated in parallel by separate uni-modal source estimations when $\alpha$ is set to 0 and 1, respectively (with $\lambda _j=1$); afterwards, the joint source estimation tunes both $\alpha$ and $\lambda _j$. This optimization algorithm is shown in Algorithms 1 and 2 for better clarification. The performances of estimations for different parameters are then illustrated in Fig. 2. The detected values of $\lambda _e^*$ and $\lambda _b^*$ respectively equal to $3.4\times 10^{-5}$ and $10^{-3}$, show the different order of magnitude of the two data-fidelity terms, and confirms our claim on the need of calibration parameters. Moreover, the value of $\alpha ^*=0.6$ maximizes the performance which shows that the estimation is almost equally balanced between the two EEG and BOLD data, and that the uni-modal parameters are well calibrated. Therefore, $\alpha$ can always be set to 0.5 and $\lambda$ is the only parameter that remains to be estimated in such method. This calibration process can thus reduce the number of parameters and simplifies the procedure.

Fig. 2.

(a) selection of $\lambda _e^*$ with EEG only performance; (b) selection of $\lambda _b^*$ with fMRI only performance; (c) selection of $\lambda _j^*$ and $\alpha ^*$ with EEG and fMRI jointly.

Fig. 3.

Box plots over eight subject of the power score (8) for activations of the union of the PMA and M1 regions in both hemispheres, while doing a right-hand motor task.

Fig. 4.

Estimated location of sources. Blue zones: active sources shown according to their power score (8). Red circles: Unexpected detections of activation.

Discussions. Here, some results of the joint method are depicted based on the estimated parameters. Figure 3 shows the detected power score, defined in Eq. (8), of the estimation in the motor cortex. Since the task brings more attention toward the activation of PMA and M1, the activities of the left and right PMA and M1 for all subjects are compared. A higher difference between the power of the left and right regions is seen in the joint approach ($p=0.02$ for a right-tailed paired t-test), while this difference is not significant in the uni-modal ones ($p\,=\,0.1$ for BOLD-only and $p=0.4$ for EEG-only).

The spatial result of a subject is illustrated in Fig. 4. In fMRI-based method, increased activity is seen in the left motor cortex. The activities in unexpected regions can be either due to vascular activities or the errors of the BOLD linear model [17]. However, a precise localization of the active areas (SMA and M1) can be seen; while EEG detected activities in SMA but not the left motor cortex. The unexpected activities can be due the fact that every instantaneous neuronal activity is recorded in EEG because of its high temporal precision. The multi-modal approach is capable of avoiding the false positives of the uni-modal approaches. Besides, it has detected high activities in both SMA and the left motor cortex.

Fig. 5.

Temporal activity of a source detected in SMA by uni- vs multi-modal approach.

Finally, the temporal signal of a detected source is shown (Fig. 5). This result cannot be quantitatively validated. Analysis can be on spectral content of EEG. Here, the temporal resolution is only compared, and further analysis is left for future studies. Here, the detected source has low temporal resolution in fMRI-only approach comparing the EEG-only one. The parts detected as zero are caused by the sparsity constraint. The temporal resolution of the estimation detected by the joint approach is as high as the temporal resolution of EEG.

5 Conclusions and Perspectives

This paper showed that a symmetric multi-modal approach based on EEG and fMRI, can improve the spatio-temporal resolution of the estimation of neuronal sources. Linear models of EEG and fMRI are used, and the sparsity constraint is applied on the sources as a logical hypothesis which helps finding solutions to the undetermined EEG and BOLD inverse problems. This model is proved to be successful by the results on real data of 8 subjects doing a right hand motor task. Moreover, the detected values of the model parameters validated the efficiency of the joint estimation when both data are symmetrically and almost equally involved in the model. This work supports our assumption that a simple joint linear model can work on real data, so that in future studies more complex models can be taken into consideration.