Keywords

1 Introduction

To infer the activated brain sources from the recorded EEG data is called inverse problem. Precise localization of neuronal firing pattern inside the brain can offer an insightful understanding of how the brain is functioning under certain cognitive and motion tasks. We also argue that source reconstruction or solving the inverse problem is the first and primary step for connectivity analysis of the brain, and precise inference of time course of brain sources is required in order to build the brain connectivity network. The latter step is to analyze the brain network using complex networks [4, 14, 17] characteristics measurement, as we saw a shift in neuroscience community from traditional “segregation” perspective to “integration” perspective where the functional and effective connectivity between different regions of brains are intensively studied [8, 13] in recent years.

In order to solve the ill-posed inverse problem, different priors or assumptions have to be imposed to obtain a unique solution. The most traditionally used priors are based on minimum power, leading to what is known as the minimum norm estimate (MNE) inverse solver [6], or minimum magnitude, termed as minimum current estimate (MCE) [16], leading to a least absolute shrinkage and selection operator (LASSO) formulation. Other assumptions or priors are presented with different inverse algorithms, such as standardized low-resolution brain electromagnetic tomography sLORETA [15], which enforces spatial smoothness of the source located on neighboring voxels; Bernoulli-Laplace priors, which introduced \(\ell _0\,+\,\ell _1\) norm in a Bayesian framework [2]; Mixed Norm Estimates (MxNE), which imposes sparsity over space and smoothness over time using \(\ell _{1,2}\)-norm regularization [3]; graph Fractional-Order Total Variation (gFOTV) [11], which impose sparsity of the spatial fractional derivatives so that it locates source peaks by providing the freedom of choosing smoothness order.

As summarized above, numerous algorithms that were based on different source configuration assumptions or prior knowledge were presented to solve the inverse problem. Traditional algorithms solve the EEG inverse problem independently for different brain states without leveraging the label information, that will make it hard to compare the reconstructed sources under different brain states due to its low SNR (Signal-to-Noise Ratio) of the EEG signal. To the best of our knowledge, few researchers come up with a model that can integrate EEG inverse problem with label information (e.g. happiness, sadness, and surprise) to find task related discriminative sources except our previous work [12]. To explicitly extract factual sources and eliminate the spurious ones, we proposed the graph regularized version of discriminative source reconstruction that has the capability of promoting intra-class consistency, and we tested it on synthetic data and illustrated its effectiveness in discovering task related sources.

2 The Inverse Problem

Under the quasi-static approximation of Maxwell’s equations, the EEG signal measurements X can be described as the following linear function of current sources S:

$$\begin{aligned} X= LS + E, \end{aligned}$$
(1)

where \(X \in {\mathbb {R}^{{N_c} \times {N_t}}} \) is the EEG data measured at a set of \(N_c\) electrodes for \(N_t\) time points, \( L \in {\mathbb {R}^{{N_c} \times {N_d}}}\) is a wide matrix called lead field matrix that maps the brain source signal to sensors on the scalp, each column of L represents the activation pattern of a particular brain source to the EEG electrodes, \(S \in {\mathbb {R}^{{N_d} \times {N_t}}}\) represents the corresponding neural electric field in \(N_d\) source locations for \(N_t\) time points. \(E \in {\mathbb {R}^{{N_c} \times {N_t}}} \) is additive noise. An estimate of S can be found by minimizing the following cost function, which is composed of a data fidelity term and a regularization term:

$$\begin{aligned} \arg \mathop {\min }\limits _S \left\| {X - LS} \right\| _F^2 + \lambda \varTheta (S), \end{aligned}$$
(2)

where \({\left\| \cdot \right\| _F}\) is the Frobenius Norm. The regularization term \(\varTheta (S)\) can be used to guarantee smooth source configurations temporally or spatially and enforces neurophysiologically plausible solutions or to guarantee sparsity in source solution. For example, to restrict the total number of activated voxels to be less than or equal to k, the constraint \({\left\| {{s_i}} \right\| _0} \leqslant k\) can be used. Even though \(\ell _0\)-norm is the best intuitive formulation to restrict number of activated sources, it’s a common practice to use approximated norm such as \(\ell _1\) to avoid the problem being NP-hard when solving EEG inverse problem. For the ith time point, the \(\ell _1\) regularized formulation is given below:

$$\begin{aligned} \left. {\left\langle {s_i} \right. } \right\rangle = \arg \mathop {\min }\limits _{s_i} \left\| {{x_i} - L{s_i}} \right\| _2^2 + \gamma {\left\| {{s_i}} \right\| _1}. \end{aligned}$$
(3)

Given the EEG recordings at a time point, which is denoted as ith column \(x_i\) of X matrix, we want to represent the signal with minimum error by trying to find the best linear representation from activation patterns (atoms) in the over-complete dictionary L [12]. The solution \(s_i\) is the sparse coding for the \(x_i\) in the dictionary L, the non-zero entries in \(s_i\) corresponding to a column in the dictionary matrix L represent the activated regions inside the brain [12].

3 Proposed Framework

3.1 Graph Regularized Discriminative Source Imaging

Due to the fact that EEG signal is non-stationary and typically the SNR is very low, it’s important to get consistent inverse solution under the same brain status and eliminate the spurious sources that are usually not consistent within the same class. Inspired by the successful applications of graph regularization in computer vision community [1, 5], the proposed model of retrieving task related discriminative source is presented, which is termed as Graph Regularized Discriminative Source Imaging (GRDSI), and comprises the data fidelity term and label guided graph regularization term:

$$\begin{aligned} {\begin{matrix} \left. {\left\langle {S} \right. } \right\rangle = \arg \mathop {\min }\limits _{S} \left\| {X - LS} \right\| _F^2 + \gamma {\left\| S \right\| _{1,1}} + \frac{\beta }{2}\sum \limits _{i,j = 1}^N {\left\| {{s_i} - {s_j}} \right\| _2^2{M_{ij}}} , \end{matrix}} \end{aligned}$$
(4)

where the first term is the fidelity term, the second term is the cost of sparse coding, \({\left\| \cdot \right\| _{1,1}}\) is the \(\ell _1\) norm notation for a matrix, equal to the sum of the absolute value of all elements in a matrix. The third term is the graph regularization term that requires all the sparse coder within the same category remains similar pattern while making the sparse representation for different class distinct from each other. The definition of M matrix can be written as:

$${M_{ij}} = \left\{ {\begin{array}{*{20}{l}} {1,\;{\text {if (}}{s_i}{\text {,}}{{\text {s}}_j}{\text {) belong to the same class}}{} \quad } \\ {0,\;{\text {if (}}{s_i}{\text {,}}{{\text {s}}_j}{\text {) belong to different classes}}} \\ \end{array}} \right. $$

The goal of this formulation is to find discriminative sources while maintaining the robustness of in-class reconstructed sources.

Remarks on design of M matrix

When \((s_i,s_j)\) belong to the same class, design the value of \(M_{ij}\) to be positive will add penalty when the intrinsic geometry \((s_i,s_j)\) is different, thus promoting intra-class consistency of the source and reduce the spurious sources estimated at each time point.

Define D as a diagonal matrix whose entries are column or row sums of the symmetric matrix M, \(D_{ii} = \sum \nolimits _j {{M_{ij}}} \), define \(G = D - M\), G is called graph Laplacian [1], The third term of Eq. 4 can be further derived in the following way:

$$\begin{aligned} \sum \limits _{i,j = 1}^N \left\| {{s_i} - {s_j}} \right\| _2^2{M_{ij}} = \sum \limits _{i,j = 1}^N ({s_i}^Ts_i+{s_j}^Ts_j -2{s_i}^Ts_j )M_{ij} = 2tr(SGS^T) \end{aligned}$$
(5)

As a result, Eq. 4 is further derived as:

$$\begin{aligned} \left. {\left\langle {S} \right. } \right\rangle = \arg \mathop {\min }\limits _{S} \left\| {X - LS} \right\| _F^2 + \gamma {\left\| S \right\| _{1,1}} + \beta (Tr(SGS^T)) \end{aligned}$$
(6)

Equation 6 can be efficiently solved using feature-sign search algorithm [1, 10].

3.2 Common Sources Decomposition with Voting Orthogonal Matching Pursue (VOMP)

Under the assumption of strong common spontaneous source activation pattern, the contribution of discriminative sources to the EEG recorded data is relatively small, making the solution space for different classes highly correlated and difficult to find discriminative sources. As a result, the convex hull spanned by all the source configuration is limited to a tiny portion of the space [18]. In order to address that, we use the idea of “cross-and-bouquet” model [18] and come up with a useful step that is to decompose of X to find the common sources shared by different classes. The Voting Orthogonal Matching Pursue (VOMP) is proposed and described in Algorithm 1. The aim is to recover the common sources across all classes. The core part of VOMP is Orthogonal Matching Pursue (OMP) which is a very efficient algorithm. After the decomposition of common source, its contribution to the EEG data X is also removed. The new EEG data after removal of the common source is written as \(X_{new} = X - LS_c\).

Based on the discussion above, the proposed framework to solve Problem 6 is summarized in Algorithm 2 and illustrated in Fig. 1.

figure a
Fig. 1.
figure 1

Procedures of our framework: After gathering labeled EEG recorded data, the brain model is constructed using finite element method (BEM) based on MRI images, the VOMP algorithm is used to decompose the primary common source starting with a high minimum voting percentage, and then solve it using feature-sign search algorithm, the last step is to map discriminative sources to the cortex.

Full size image
figure b

4 Numerical Results

We used a recently developed realistic head model called ICBM-NY or “New York Head” [9]. The dimension of lead field matrix we are using is \(108 \times 2004\), representing 108 channels and 2004 voxels. We also assume that source orientation is perpendicular to the cortex surface. In each simulation, noises originate from sensor level and cortex voxel level both contributed to the recorded EEG data. The SNR is calculated as \(SNR = 20\log _{10} \frac{{\left\| S \right\| _2}}{{\left\| N \right\| _2}}\). We show the effectiveness of the graph regularization term in reconstructing the discriminative sources by comparing it with the other eight benchmark algorithms, including ElasticNet, Homotopy, DALM, PDIPA, L1LS, FISTA, sLORETA, MNE. The former 6 algorithms are compared in image reconstruction applications and can be referred to Reference [19] for details.

We designed the spontaneous common sources with a magnitude of 0.8 with standard deviation to be 0.1 and task related discriminative source with a magnitude of 0.2 with a standard deviation of 0.05 located in different Region Of Interest (ROI)s from the common sources. The ROI we used here are defined in Reference [7]. We sampled 200 time points for each class and did the experiment 5 times to get the average accuracy of the reconstructed source. For the GRDSI parameter, we set \(\beta \) to be 0.05 and \(\alpha \) to be 0.06; The noise matrix is designed to affect the EEG recording together with the true source signal. For each time point, 3 random voxels are corrupted randomly with the average value being 0.2, 0.4, 0.6 and variance being 0.05 based on different SNR design. All computations were conducted on a 64–bit Linux workstation with 3.00 GHz i7-5960x CPU and memory of 64 GB.

The reconstruction performance of the proposed method as well as the benchmark methods based on 150 experiments are summarized in Table 1. All of the values in Table 1, except the Time column (in seconds) represents distance in (mm) from ground true source to the reconstructed source calculated from the shortest path along cortex surfaces. PSE represents primary source error, which is the distance of reconstructed primary source to the ground truth primary source. PSE measures the capability of each algorithm to reconstruct the common sources. When the reconstructed location is in a different hemisphere from the ground truth, there is no path connecting those two voxels, so we mark the distance to be 250 mm. EC1 represents error for class 1, which is the distance of the reconstructed discriminative source to the ground truth. EC2 and EC3 are similarly defined.

Fig. 2.
figure 2

Ground truth for all 3 classes

Full size image
Table 1. Reconstruction accuracy summary
Full size table
Fig. 3.
figure 3

MNE solution: The above row is the MNE solution for class 1; Class 2 and class 3 is illustrated in the middle and bottom row. The solution MNE gives is not sparse, with too many spurious sources of small magnitude.

Full size image
Fig. 4.
figure 4

sLORETA inverse solution: sLORETA can successfully reconstruct the primary source, however the secondary source is not successfully reconstructed. Compared to the solution of MNE, sLORETA can suppress spurious sources with small magnitude.

Full size image
Fig. 5.
figure 5

GRDSI reconstructed source: The reconstruction solutions for 3 classes are given in each row. The discriminative source can be successfully reconstructed.

Full size image

To illustrate the effect of the proposed framework, the ground truth of the activated pattern is given in Fig. 2, with the reconstructed source by MNE, sLORETA and our method given in Figs. 3, 4 and 5. We can see from Table 1 and the Figs. 3, 4 and 5 that when the SNR is large, all the algorithms performs well in reconstructing primary source, as for the discriminative sources for different classes, our method can achieve almost perfect reconstruction. All other algorithms’ performances are also acceptable when SNR is large, except for sLORETA, MNE and ElasticNet. When we increase the noise, all of the algorithms can still achieve high accuracy in finding the primary source. For the discriminative source, our algorithm performs much better. We also validated that, to solve a pure \(\ell _1\) EEG inverse problem, the Homotopy algorithm performs better in most cases than other algorithms in the EEG inverse problem, which is in line with Reference [19].

5 Conclusion

In this paper, we proposed to use label information to retrieve discriminative sources corresponding to different brain status. A graph regularized EEG inverse formulation called GRDSI that implicitly uses the label information was presented that can boost the intra-class consistency and eliminate spurious sources. We bring up the idea of “cross-and-bouquet” in the inverse problem and present an efficient algorithm to address the high coherence problem of the reconstructed signals given high background spontaneous source signal. An efficient algorithm called feature-sign search algorithm is used to solve the GRDSI problem. We illustrated the superior of our algorithm in retrieving discriminative sources while traditional algorithms failed given certain level of noises.