The relationship between the scalp potential measured by the electrodes and the brain source distribution can be expressed as follows:

where t represents the time, vector x(t) ∈ R^{n × 1} represents the EEG or MEG signal measured by n electrodes, matrix H ∈ R^{n × m} represents the lead field, vector s(t) ∈ R^{m × 1} represents the source signal generated by m brain sources, and vector ε(t) ∈ R^{n × 1} represents the additive noise from observation.

The forward model models the linear mapping between scalp potential measured by the electrodes and the brain source signal (Birot et al., 2014). The solution to the forward problem relies on the establishment of the head model, which is determined by the geometry and corresponding electrical conductivity of different head tissues such as brain, skull, scalp, etc. (Acar and Makeig, 2010). In the early days, the mainly used head models are the spherical model and the ellipsoid model (Gutiérrez et al., 2005). With the development of brain imaging technology, real head models, which can be calculated by the boundary element method (BEM), the finite elements method (FEM), and the finite difference method (FDM) are increasingly used (Akalin-Acar and Gençer, 2004). Once the head model is established, the lead field matrix can be determined.

The connectivity relationship between m sources can be represented by an undirected graph, which can be defined as follows:

where V ∈ R^{m × 1} is a set of m nodes, A ∈ R^{m × m} is the corresponding adjacency matrix. If there is no edge connecting nodes i and j, then a_ij = 0; otherwise, a_ij > 0, and its value represents the weight of the edge between the two nodes. In addition, since G is an undirected graph, then a_ij = a_ji, which means the adjacency matrix A is symmetric. In the EEG source space, all the potential source locations are represented by the nodes defined on a 3D mesh as is illustrated in Figure 1. When the source locations i and j are neighbors on the 3D mesh, then we set a_ij > 0.

The graph signal is defined on the set of graph nodes V, which is represented by a vector, and each element represents the signal value at the corresponding node. The brain source signal s = [s₁, s₂, …, s_m]^T, in which each element s_i represents the signal value of the i-th source voxel, is defined on the nodes of graph G. The traditional Fourier transform calculates the projection of a function f(t) on the basis function e^−iwt, and the projected value of a time series signal using Fourier transform F(w) represents its magnitude at the basis of a specific frequency. Different from the traditional Fourier transform defined in the temporal domain, for signals defined on a graph, the eigenvectors of the Laplacian matrix L of the graph can be used as the basis vectors of the GFT, where the Laplacian matrix L can be calculated as follows:

where L ∈ R^{m × m}, and D ∈ R^{m × m} is a diagonal matrix called the degree matrix, in which the diagonal elements satisfy dii=∑jmaij, that is, the sum of elements in the i-th row of A. Since A and L is real and symmetric, therefore, L can be decomposed as follows:

where Λ ∈ R^{m × m} is a diagonal matrix and the diagonal elements λ_i, (i = 1, 2, …, m) are the eigenvalues of L and satisfy λ₁ ≤ λ₂ ≤ … ≤ λ_m, U ∈ R^{m × m} is the eigenvector matrix, which is also an orthogonal matrix satisfying UU^T = I, each column in U is an eigenvector of L and corresponds to the eigenvalue in Λ. With the eigenvectors of the Laplacian matrix L as the Fourier basis vectors, each eigenvector can be regarded as a graph basis with a certain frequency, and this frequency corresponds to the eigenvalue. The smaller the eigenvalue, the lower the frequency of the corresponding eigenvector, which is manifested as a small difference between the signals of adjacent nodes on the graph; on the contrary, a larger variation among neighboring signals. The value of a graph Fourier coefficient can measure the amplitude of the graph signal at different frequencies. With the eigenvectors as the Fourier basis vectors, the GFT of a given graph signal s can be defined as follows:

where vector s∼∈Rm×1 is the graph Fourier coefficient. Further, the inverse graph Fourier transform (IGFT) of s can be defined as:

The above two formulas show that a graph signal can be decomposed into components with different frequencies through the GFT, and can also be recovered through the IGFT.

To characterize the graph frequency, we introduce the following definition:

Definition 1: Graph Frequency (GF): GF, denoted as f_G, is a function of u_i which represents the total number of sign flips of u_i between any two connected nodes on G, it is defined as follows:

where Ω(j) represents all neighbors of node j, and I(⋅) is an indicator function to check whether the values of u_i at nodes j and p have a sign flip. The number of sign flips is analogous to counting the number of zero crossings of the basis signal within a given window for a time series data. We constructed the Laplacian matrices within first-order neighbors and second-order neighbors, and the associated GF spectrum is shown in Figure 2. It can be seen that the GF value of the eigenvector increases as the eigenvalue increases.

Similar to the counterpart in the time domain, the spatial frequency basis matrix U can be similarly decomposed into different spatial frequency bands, such as U = [U_low, U_medium, U_high]. The graph signal s can be projected into a subspace of U. For example, s∼=UlowTs is the projection of s into a space spanned by low frequency eigenvectors. In our work, we use the spatially low frequency components as a filter to reconstruct the focally extended sources.

Bidirectional long-short term memory neural network is an extension of the traditional RNN (Schuster and Paliwal, 1997). For the time series, it is recognized that RNN can effectively estimate the information at the future moment based on the previous states. However, for the time series with a long sequence of states, the estimation performance of RNN will be greatly discounted, because the future information in a long time series usually depends on the information from distant history moments, which is the long-term dependence. However, the superior structure of the LSTM unit equips the network with the ability to solve long-term dependence. The RNN with LSTM units can filter the information by a unique structure called “gate” and store the valid information by the so-called “memory cell.” The elements in a gate vector have values in the interval [0,1]. When preceding time series information arrives at the gate, it will be multiplied with the gate vector element-wise, if the element value in the gate vector is 1, then the timing information multiplied with it will be retained, and if the element value is 0, the information will be discarded after a multiplication with 0. In this way, the filtering of information propagated in the LSTM unit is achieved. The valid information obtained after filtering is then stored in the memory cell and passed on to the next moment to prevent being lost over time, thus effectively addressing the long-term dependence existing in traditional RNNs.

The structure of a standard LSTM unit is shown in Figure 3A, where x_t, h_t−1, and c_t−1, respectively, represent the input sample at the current moment, the unit output and the memory state at the previous moment. The information contained in x_t and h_t−1 is first activated by σ(⋅) function and got the forget gate f_t, the input gate i_t, and the output gate o_t. At the same time, x_t and h_t−1 are also activated by tanh(⋅) function and got a temporary state c∼t. On the one hand, the information passed from the previous moment which contained in c_t−1 is filtered by the forget gate f_t. On the other hand, the newly input information contained in c∼t is filtered by the input gate i_t. Then, the valid information retained by the above two filtering processes is integrated together as a new memory state c_t. This newly updated memory state is passed along time to the next moment, and simultaneously, it is also filtered by the output gate o_t. Finally, the new output of the LSTM unit h_t is obtained.

This propagation process can be formulated as follows:

where W_f, W_i, W_c, W_o are weight matrices; b_f, b_i, b_c, b_o are bias vectors; the symbol * stands for the element-wise multiplication.

The hidden layer of the BiLSTM neural network is composed of two layers of LSTM units that are reversely connected, and its structure is shown in Figure 3B. In a BiLSTM layer, the time series performs both forward propagation and backward propagation. Therefore, both the information at the previous moments and that at the future moments can be fully utilized.

The output of the BiLSTM neural network can be calculated as follows:

where y_t is the final output, ht′ is the unit output of the backward propagation, W_s and b_s are the weight matrix and the bias vector of the output layer, respectively, and the symbol ⊕ stands for the vector concatenation.

To render source extents activation, the adjacent sources along with a central source are activated at the same time. The signal strength of the adjacent sources is set to be lower than that of the central region. All the 2,052 potential source locations are chosen as the central source in turn to generate the scalp EEG data. In the first experiment, we test the proposed algorithm on the simulated EEG data and the true source activation pattern with one source extents activated. Apply the training and validation dataset to train and validate the proposed GFT-BiLSTM model, and then test it on the testing dataset. The performance of the GFT-BiLSTM is compared with that of dSPM, MNE, and sLORETA. The evaluation metrics of each model are shown in Table 1 and Figure 5. The comparison between the ground truth source and the reconstructed sources by different algorithms is shown in Figure 6.

From Table 1 and Figure 5, it can be seen that the proposed GFT-BiLSTM shows the better performance when compared to other methods. For different SNR levels, the AUC corresponding to GFT-BiLSTM is the highest while the LE is the lowest. The brain source distribution estimated by the proposed GFT-BiLSTM is closer to the ground truth as illustrated in Figure 6. The numerical result demonstrates the superiority of the GFT-BiLSTM when applied to solve the ESI inverse problem.

In order to study the performance of the proposed GFT-BiLSTM when there are multiple activated sources, we randomly select 2 out of 2,052 brain source locations to be activated and the first level neighboring sources of these two source locations are also activated with a lower signal magnitude. Apply this simulated dataset to train and validate the GFT-BiLSTM, and then test it on the testing dataset. The performance of the GFT-BiLSTM is compared with that of dSPM, MNE, and sLORETA. The evaluation metrics of each model are shown in Table 2 and Figure 7. The comparison between the real source and the estimated source is shown in Figure 8.

From Table 2 and Figure 7, the proposed GFT-BiLSTM demonstrates better performance when compared with other methods. In addition, with the increased number of activated source, the estimation performance of all methods except the GFT-BiLSTM deteriorates significantly. This is demonstrated as an increase in LE and a decreased AUC value. There is also a situation for other method in which the localization of the central region is accurate while the localization of its adjacent regions is slightly deviated. The reason is that as the number of activated regions increases, the distribution of the sources is no longer concentrated, and it is more challenging to accurately estimate the locations of all active regions. The reduction in performance for AUC is much more pronounced for the benchmark algorithms. In contrast, the performance of the GFT-BiLSTM is more stable and robust when it comes to multiple activated sources.

In order to further evaluate the performance of the proposed GFT-BiLSTM, we applied it on the public epilepsy EEG dataset from the Brainstorm tutorial datasets (Tadel et al., 2011). This dataset was recorded from a patient who suffered from focal epilepsy. The patient underwent invasive EEG to identify the epileptogenic area then underwent a left frontal tailored resection and was seizure-free during a 5-year follow-up period. We followed the Brainstorm tutorial to obtain the head model, and the lead field matrix. Then we calculated the average spikes (as shown in Figure 9) of the provided EEG measurements with 29 channels. Apply the averaged EEG data for brain source localization, and the estimated sources at peak (0 ms) from different methods are shown in Figure 10, as compared to other methods including dSPM, MNE, sLORETA.

It can be seen from Figure 10 that the proposed GFT-BiLSTM provides a good reconstruction of the epileptogenic zone which was validated by the follow-up survey after resection on the left frontal region. The source area estimated by dSPM and sLORETA spans a wide range cortical areas and includes part of the right frontal lobe which is not related to the epilepsy lesion. In contrast, the source location estimated by the MNE method and the GFT-BiLSTM proposed in this paper is more accurate. However, compared between the two methods, the range of sources estimated by the GFT-BiLSTM is smaller, and the source estimated by GFT-BiLSTM shows better continuity of the spatial signal, due to the benefit of using GFT.