Rejuvenating classical brain electrophysiology source localization methods with spatial graph Fourier filters for source extents estimation

Abstract

EEG/MEG source imaging (ESI) aims to find the underlying brain sources to explain the observed EEG or MEG measurement. Multiple classical approaches have been proposed to solve the ESI problem based on different neurophysiological assumptions. To support clinical decision-making, it is important to estimate not only the exact location of the source signal but also the extended source activation regions. Existing methods may render over-diffuse or sparse solutions, which limit the source extent estimation accuracy. In this work, we leverage the graph structures defined in the 3D mesh of the brain and the spatial graph Fourier transform (GFT) to decompose the spatial graph structure into sub-spaces of low-, medium-, and high-frequency basis. We propose to use the low-frequency basis of spatial graph filters to approximate the extended areas of brain activation and embed the GFT into the classical ESI methods. We validated the classical source localization methods with the corresponding improved version using GFT in both synthetic data and real data. We found the proposed method can effectively reconstruct focal source patterns and significantly improve the performance compared to the classical algorithms.

1 Introduction

EEG/MEG are non-invasive measurement modalities with high temporal resolution up to 1 ms. EEG/MEG signals can be collected noninvasively through electrodes or sensors on the scalp. Importantly, EEG/MEG are the direct measurement methods to detect the instantaneous electrophysiological activities of the brain [1]. EEG and MEG are recognized as powerful tools for capturing real-time brain functions by measuring neuronal processes with broad clinical and neuroscience applications [2]. The EEG/MEG source localization aims to solve an inverse problem to reconstruct the underlying electrophysiological brain activities given the measured EEG/MEG signal and the brain forward model [3]. The EEG/MEG source localization is also referred to as Electrophysiological Source Imaging (ESI) [4]. However, the number of EEG/MEG channels is far less than that of the brain sources, which makes the ESI an ill-posed problem.

In the past decades, numerous algorithms have been developed with different assumptions on the configuration of the source signal. One seminal work is minimum norm estimate (MNE) where $\ell _2$ norm is used as a regularization [5], which is to explain the observed signal using a potential solution with the minimum energy. Different variants of the MNE algorithm include dynamic statistical parametric mapping (dSPM) [6] and standardized low-resolution brain electromagnetic tomography (sLORETA) [7]. The $\ell _2$-norm-based methods tend to render spatially diffuse source estimation. To promote sparse solutions, Uutela et al. [8] introduced the $\ell _1$-norm, known as minimum current estimate (MCE). Rao and Kreutz-Delgado proposed an affine scaling method [9] for a sparse ESI solution. The focal underdetermined system solution (FOCUSS) proposed by Gorodnitsky et al. encourages a sparse solution by introducing the $\ell _{p}$-norm regularization [10]. Besides, Bore et al. proposed to use the $\ell _{p}$-norm regularization ($p < 1$) on the source signal and the $\ell _1$ norm on the data fitting error term to achieve sparsity [11]. Babadi et al. demonstrated that sparsely distributed solutions to event-related stimuli could be found using a greedy subspace-pursuit algorithm [12]. Wipf et al. proposed a unified Bayesian learning method that can automatically calculate the hyperparameters for the inverse problems under an empirical Bayesian framework, and the sparsity of the solution is also guaranteed [13]. To mitigate the high level of noise when estimating source signal, Hashemi et al. applied a hierarchical Bayesian (type-II maximum likelihood) model for joint estimation of latent variables for both source and noise [14]. An innovative robust empirical Bayesian framework proposed by Ghosh et al. can estimate the low-rank noise covariance and sparse brain source activity at the same time [15]. Wan et al. reformulate the ESI problem into a graph search problem by exploiting the graph neighborhood information in the brain source space and the optimal solution can be theoretically guaranteed under high noise circumstances [16].

In addition to using optimization approaches, deep learning methods have become increasingly popular in recent years. Hecker et al. came up with a CNN model with a special design of the input matrices of EEG and MEG signal, termed as ConvDip, for source localization which achieved good performance [17]. Jiao et al. proposed to use of a graph Fourier transform based bidirectional long-short-term memory, termed as GFT-BiLSTM, which can reduce the output dimensions and can construct an extended brain region activation. The GFT-BiLSTM method achieves good performance localizing seizure onset zone [18]. To take advantage of both the high spatial resolution of fMRI and the high temporal resolution of EEG, Liu et al. proposed a hierarchical deep transcoding model for fusing simultaneous EEG-fMRI data to estimate the brain source activity with high spatiotemporal resolution [19, 20]. Another multi-modal deep learning model that fuses both MEG and EEG information is shown to have high accuracy compared to using a single modality of EEG or MEG in the context of deep learning approach [21].

As the brain sources are not activated discretely due to its volume conductivity property, a compact and extended area of source estimation is preferred [4, 22, 23], and it has been used for multiple applications, such as somatosensory cortical mapping [24], and epileptogenic zone in focal epilepsy patients [25, 26]. To estimate an extended area of source activation, we propose to use the GFT technique, which has been shown with an improved accuracy in reconstructing extended source activations [18, 27]. GFT plays a pivotal role in understanding the spectral characteristics of the signals on a graph, analogous to how the Fourier Transform reveals the frequency content of a signal in the time domain [28]. By decomposing graph data into its spectral components, it provides a unique perspective on the underlying structure and connectivity/continuity patterns within the graph. Numerous interesting applications can be found in various domains. For example, in their seminal work, Defferrard et al. proposed to use GFT with deep learning and came up with ChebNet, which highlighted the use of GFT in neural networks for graph-structured data [29]. GFT has been used to characterize the spectrum properties of the brain connectivity networks [30], Brahim and Farrugia showed that using GFT-based structural connectivity and functional connectivity analysis can provide a more accurate classification for patients with Autism Spectrum Disorders [31].

In this work, we propose to improve the classical ESI methods using GFT and its low-frequency components as source space subspaces, and we highlight the importance of using GFT for an extended area of source activation before applying the classical source imaging methods. A preliminary version has been published in the Brain Informatics Proceedings [32].

2 Method

In this section, we start with the introduction of the ESI inverse problem, which is followed by the presentation of the GFT and the improved version of classical methods using GFT.

2.1 EEG/MEG source imaging

The source imaging forward problem can be described in the formula as $Y = KS + E$, where $Y\in \mathbb {R}^{C\times T}$ is the EEG or MEG measurements, C is the number of EEG or MEG channels, T represents the time sequence length, $K\in \mathbb {R}^{C\times N}$ is the leadfield matrix which is a linear mapping from the brain sources to the EEG/MEG sensors, N is the number of brain source, $S\in \mathbb {R}^{N\times T}$ represents the brain source signal in N source locations for all the T time points, and E is the measurement noise. The inverse problem is to estimate S given Y and K. Since the source dimension N is much larger than the number of electrodes C, making the ESI an ill-posed inverse problem, it is challenging to obtain a unique solution. Various regularization terms were designed based on the prior assumptions of the spatial and temporal structure of the source signal to find a unique solution. The inverse problem of ESI can be formulated as below:

$$\begin{aligned} S = {\mathop {\textrm{argmin}}\limits _S}\, \frac{1}{2} \left\| {Y - KS} \right\| _F^2 + \lambda R (S), \end{aligned}$$

(1)

where ${\Vert \cdot \Vert _F}$ is the Frobenius norm, and S can be obtained by solving the minimizing problem. The first term in Eq. (1) is the data fitting term trying to explain the recorded EEG/MEG measurements. The second term is called the regularization term, which is imposed to find a unique solution using sparsity or other neurophysiology-inspired regularizations. For example, if R(S) is a $\ell _2$ norm, the problem is called minimum norm estimate (MNE).

2.2 Graph Fourier transform in brain source space

Consider an undirected graph $G = \left\{ \mathcal {V}, A\right\}$ generated from the 3D mesh of brain cortex, where $\mathcal {V} = \left\{ v_1, v_2, \ldots , v_N \right\}$ is the set of N nodes, A is the weighted adjacent matrix with entries given by the edge weights $a_{ij}$ that represents the connection strength between node i and node j. The graph Laplacian matrix is defined as $L = D - A$, where D is the degree matrix with $D_{ii} = \sum _{j\ne i}A_{ij}$. Since L is a positive semi-definite matrix, its eigenvalues are all greater or equal to 0, and the associated eigenvectors $U = [u_1, u_2, \ldots , u_N]$, $U \in \mathbb {R}^{N \times N}$ can be regarded as the basis vectors of GFT where any signal in the graph can be approximated as the linear combinations of basis. Thus, the GFT for a signal S can be defined as $\tilde{S}= U^{T}S$, whereas the inverse GFT is given as $S = U\tilde{S}$. Here we define normalized graph frequency (NGF) [27] as

$$\begin{aligned} f_G(u_i) = \frac{f_s(u_i)}{Tr(L)}, \end{aligned}$$

(2)

where Tr(L)is the trace of L, and $f_s(u_i)$ is defined as

$$\begin{aligned} f_s(u_i) = \sum _{m=1}^{N}\sum _{n\in \mathcal {N}(m)}\mathbb {I}(u_i(m)u_i(n)<0)/2, \end{aligned}$$

(3)

where $\mathcal {N}(m)$ represents all neighbors of node m, and $\mathbb {I}(\cdot )$ is the indicator function which equals 1 if the values of $u_i$ on node m and n have different signs and 0 otherwise. The number of sign flips at time t indicates how many zero crossings of a signal are within a bounded region at t.

We calculated the NGF in the whole time series within first-order neighbors, second-order neighbors, and third-order neighbors, respectively. The spectrogram, which is illustrated in Fig. 1, reveals that the NGF is positively correlated with the order of the eigenvalue of L. Thus we can further separate U into low, medium, and high-frequency components according to NGF values, and reformat it as $U = [U_{L}, U_{M}, U_{H}]$. The source activation patterns corresponding to eigenvectors with different NGFs are illustrated in Fig. 1.

Fig. 1

Source distributions corresponding to eigenvectors with different NGFs

2.3 Brain source imaging with spatial graph filters

The existing sparse source localization methods with $\ell _1$-norm and $\ell _{2,1}$-norm as regularizations usually provide over-sparse solutions when estimating source extents, whereas the other non-sparsity methods commonly result in multiple foci. In this work, we propose to use the low-frequency subspaces spanned by the low-frequency graph filters (eigenvectors) to approximate an extended area of brain source localization. The proposed GFT technique can improve the performance of classical source localization methods by providing the low-frequency spatial graph filters $[U_{L}] \in \mathbb {R}^{N\times P}$, thus removing the high-frequency noises to reconstruct the focally extended sources. Here P serves as the cutoff value for the set of low-frequency components in the subspace. The intuition is that the main energy of the source signal usually lies in the low-frequency components which are associated with the regions on the cortex with relatively large source extend area in a time series. Keeping the low graph frequency could promote a source extend area reconstruction and decrease the impact of the noise. Moreover, the reduced dimensional estimation in the inverse problem could further constrain the solution space and make the solution more easily solved and robust.

Specifically, we project graph signals S onto a subspace spanned by $\tilde{U}$ with low NGF values. The data fitting term can be written as $\Vert Y- K\tilde{U}\tilde{U}^T S\Vert ^2_F$. Let $\tilde{S} = \tilde{U}^T S$ and $\tilde{K} = K \tilde{U}$, then the data fitting term can be rewritten as $\Vert Y- \tilde{K}\tilde{S}\Vert ^2_F$. Solving S is equivalent to finding parameters in the subspace $\tilde{U}$. Our goal is to estimate $\tilde{S}$, and all regularizations will be added to $\tilde{S}$ instead of S.

$$\begin{aligned} \hat{S}^*_{GFT} = {\mathop {\textrm{argmin}}\limits _{\tilde{S}^*}}\, \frac{1}{2} \Vert {Y - \tilde{K}\tilde{S}^*} \Vert _F^2 + \lambda R (\tilde{S}^*), \end{aligned}$$

(4)

Finally, the source estimation $\hat{S}_{GFT}$ can be simply obtained from the inverse GFT $\tilde{U}\hat{S}^*_{GFT}$.

We delineate the solutions for 5 classical ESI methods, including MNE [33], MCE [8], MxNE (with L21 implementation) [34], dSPM [6] and sLORETA [7], and show the solution or objective function of both original form and GFT-based form as follows.

MNE MNE is proposed by Hämäläinen by imposing an $\ell _2$ regularization on the source solution [33], the MNE algorithm has a closed-loop solution, given as

$$\begin{aligned} \hat{S}&= K^T(KK^T + \lambda I)^{-1} Y, \end{aligned}$$

(5)

We define $\tilde{K}=K\tilde{U}$, then the GFT-MNE solution is

$$\begin{aligned} \hat{S}_{GFT} = \tilde{U}\tilde{K}^T(\tilde{K}\tilde{K}^T + \lambda I)^{-1} Y. \end{aligned}$$

(6)

sLORETA sLORETA is a method proposed by Pascual-Marqui [35], and the solution is also in a closed form

$$\begin{aligned} \hat{S} = {(K^T(KK^T + \lambda I)^{-1}K)^{-\frac{1}{2}}} K^T(KK^T + \lambda I)^{-1} Y, \end{aligned}$$

(7)

The GFT-sLORETA solution is given as below:

$$\begin{aligned} \hat{S}_{GFT}= \tilde{U} {(\tilde{K}^T(\tilde{K}\tilde{K}^T+ \lambda I)^{-1}\tilde{K})^{-\frac{1}{2}}} \\{} & {} \tilde{K}^T(\tilde{K}\tilde{K}^T + \lambda I)^{-1} Y. \end{aligned}$$

(8)

dSPM dSPM is proposed by Dale et al. It has a closed-form solution, and the solution is

$$\begin{aligned} \hat{S}= {(K^T(KK^T + \lambda I)^{-1}I(K^T(KK^T + \lambda I)^{-1})^T)^{-\frac{1}{2}}} \\{} & {} K^T(KK^T + \lambda I)^{-1} Y, \end{aligned}$$

(9)

the GFT-dSPM solution can be obtained as:

$$\begin{aligned} \hat{S}_{GFT}= \,& {} \tilde{U} (\tilde{K}^T(\tilde{K}\tilde{K}^T + \lambda I)^{-1} \\{} & {} I(\tilde{K}^T(\tilde{K}\tilde{K}^T + \lambda I)^{-1})^T)^{-\frac{1}{2}} \\{} & {} \tilde{K}^T(\tilde{K}\tilde{K}^T + \lambda I)^{-1} Y. \end{aligned}$$

(10)

MCE MCE is a method proposed by Uutela et al. with $\ell _1$ regularization [8]. There is no closed-loop solution. Multiple iterative solvers can solve this problem. The objective function is:

$$\begin{aligned} \hat{S} = {\mathop {\textrm{argmin}}\limits _S}\, \frac{1}{2} \left\| {Y - KS} \right\| _F^2 + \lambda \Vert S\Vert _1, \end{aligned}$$

(11)

and the GFT-MCE solution can be obtained by solving the following equation:

$$\begin{aligned} \hat{S}_{GFT}^*= \,& {} {\mathop {\textrm{argmin}}\limits _S}\, \frac{1}{2} \left\| {Y - \tilde{K}S} \right\| _F^2 + \lambda \Vert S\Vert _1 \\ \hat{S}_{GFT}= & {} \tilde{U}\hat{S}_{GFT}^*. \end{aligned}$$

(12)

MxNE MxNE used a mixed norm as a regularization term, where $\ell _2$ norm is applied to the temporal domain and $\ell _p$ ($p=\{0.5,1\}$ is applied in the spatial domain, which is proposed by Gramfort et al. [34]. There is no closed-loop solution for this regularization. In this paper, we use $\ell _{2,1}$ norm version for MxNE. The objective function is:

$$\begin{aligned} \hat{S} = {\mathop {\textrm{argmin}}\limits _S}\, \frac{1}{2} \left\| {Y - KS} \right\| _F^2 + \lambda \Vert S\Vert _{2,1}, \end{aligned}$$

(13)

and the GFT-L21 solution is given as:

$$\begin{aligned} \hat{S}_{GFT}^*= & {} {\mathop {\textrm{argmin}}\limits _S}\, \frac{1}{2} \left\| {Y - \tilde{K}S} \right\| _F^2 + \lambda \Vert S\Vert _{2,1} \\ \hat{S}_{GFT}= & {} \tilde{U}\hat{S}_{GFT}^* \end{aligned}$$

(14)

where $\Vert \cdot \Vert _{2,1}$ denotes $\ell _2$ norm on the horizontal direction and $\ell _1$ on the vertical dimension of a matrix.

3 Result

In this section, we conducted numerical experiments to validate the effectiveness of the proposed method on synthetic EEG data under different levels of neighbors (LNs) and signal Noise Ratio (SNR) settings for source localization followed by a causal time series reconstruction simulation. Lastly, we further validate it on real MEG recordings from a visual-auditory test.

3.1 Simulation experiments

To validate the proposed GFT-based methods, We first conducted experiments on synthetic data with known activation patterns.

Forward model To generate synthetic EEG data, we build a realistic head model based on T1-MRI. The brain tissue segmentation and source surface reconstruction were conducted using FreeSurfer [36]. Then a three-layer boundary element method (BEM) head was built based on these surfaces, given the tissue conductivity values. A 128-channel BioSemi EEG cap layout was used, and the EEG channels were co-registered with the head model using Brainstorm [37] and then further validated on the MNE-Python toolbox [38]. The source space contains 1026 sources in each hemisphere, with 2052 sources combined, resulting in a leadfield matrix L with a dimension of 128 by 2052.

Synthetic data generation We randomly generated 200 source locations out of 2052 locations in the source space to conduct the experiments. Furthermore, as illustrated in Fig. 2, we used three different neighborhood levels (1-, 2-, and 3-level of the neighborhood) to represent different sizes of source extents, then we activated the whole “patch" with different neighborhood levels at the same time. The activation strength of the 1-, 2-, and 3-level adjacent regions was successively set to be 80%, 60%, and 40% of the central region. The scalp EEG data was generated based on the forward model under different measurement noise configurations with different Noise Ratio (SNR) set to be 40 dB, 30 dB, 20 dB, and 10 dB. SNR is defined as the ratio of the signal power $P_{\text {signal}}$ to the noise power $P_{\text {noise}}$: $\text {SNR}=10\log (P_{\text {signal}}/P_{\text {noise}})$. In total, there were 12: 3 (levels of neighborhood) $\times$ 4 (SNRs) data sets (Y and S pairs).

Fig. 2

Brain source distributions with different levels of neighbors (LNs)

For causal time series reconstruction, we generated a simulated causal time series in the source space based on the Berlin Brain Connectivity Benchmark (BBCB) [39]. We slightly simplified it by using the autoregression with randomly generated 2 by 2 state transition matrix $\varvec{\Phi }$, with the order K to be 1 and the number of activated brain regions is set to be 2, to generate the source signal, given in Eq. 15.

$$\begin{aligned} x_t = \sum _{k=1}^K \varvec{\Phi }_k x_{t-k} + n_t, \ t=1,\ldots ,T,\ x_t\in \textbf{R}^N \end{aligned}$$

(15)

Then an independent random Gaussian noise with SNR from 10 to 40 dB was added at each time point. Lastly, a third-order Butterworth filter with zero phase delay was applied with pass bandwidth [0.1 Hz, 40 Hz].

Experimental settings We adopted MNE [5], MCE [8], $\ell _{2,1}$(MxNE) [34], dSPM [6], and sLORETA [7], as benchmark algorithms for comparison. We separately performed EEG source localization based on benchmark algorithms with and without the proposed GFT-based dimensionality reduction method. All the experiments were conducted on a Linux environment with CPU Intel(R) Xeon(R) Gold 6130 CPU @2.10 GHz and 128 GB memory.

Source localization simulation results The performance of each algorithm was quantitatively evaluated based on the following metrics:

• (1) Localization error (LE): it measures the Euclidean distance between centers of two source locations on the cortex meshes.

• (2) Area under curve (AUC): it is particularly useful to characterize the overlap of an extended source activation pattern.

Better performance for localization is expected if LE is close to 0 and AUC is close to 1. The performance comparison between the proposed methods and benchmark algorithms on LE and AUC is summarized in Table 1, and the boxplot figures, as well as the difference significance test for SNR = 30 dB, and 10 dB are given in Fig. 3. The comparison between the reconstructed source distributions with a 3-level of the neighborhood and 40 dB SNR is shown in Fig. 4.

Table 1 Performance evaluation

Fig. 3

Significance of difference through t test on AUC and LE with 3-level of the neighborhood for SNR = 30 dB (subplot A and B), and SNR = 10 dB (subplot C and D). Asterisks indicate the results of paired-sample t tests: p value < 0.5(*), p value < 0.05(**), p value < 0.005(***)

Fig. 4

Brain sources reconstruction by different ESI algorithms with the single activated area and 3-level of the neighborhood for SNR = 40 dB

From Table 1 and Figs. 3 and 4, we can find that: (1) MNE, MCE, $\ell _{2,1}$, sLORETA, and dSPM can only reconstruct the brain sources when the activated area is small and the SNR level is high, and even the evaluation metrics can be good in some cases, the reconstruction for source extend area is discrete or disconnected as shown in Fig. 4. As the source range expands and the SNR decreases, a significant increase in LE and an obvious reduction in AUC can be observed. The reconstructed source distributions are no longer concentrated. (2) By contrast, the results of the proposed GFT-based methods outperform benchmark methods in most cases after applying the spatial graph filters. They both show good stability for varied neighborhood levels and SNR settings. Particularly, the performance of the proposed GFT-based $\ell _1$ regularization family (i.e., $\ell _1$-norm and $\ell _{2,1}$-norm) exhibits better performance on reconstructing source extents without losing its advantage in sparse focal source reconstruction and outperforms other methods in most instances.

Causal time series reconstruction An order K = 1 causal time series was generated at different levels of noise. As the reconstructed time series at the source space is impacted by the strength of regularizations, to rectify this impact, we use the correlation of two sets of time series to quantify the similarity between the reconstructed time series and the ground true time series. A detailed comparison of different methods is given in Table 2, and an example of source reconstruction at 20 dB SNR is illustrated in Fig. 5.

Table 2 Pearson correlation for causal time series reconstruction

Fig. 5

Time series reconstruction on activated region at SNR = 20 dB, where the sub-figures (a)–(e) on the left column are the reconstructions of the original and proposed GFT-based forms of MCE, L21, MNE, sLORETA, and dSPM at casual source 1, whereas the sub-figures (f)–(j) on the right column are the same definition at casual source 2. The blue curves are true source signals, and the orange as well as green curves are original and proposed GFT-based forms of different methods respectively

It is worth noting that the proposed GFT-based methods are better than their original counterparts in most cases. Even at low SNR levels, the proposed methods can also perform well. Particularly, the GFT-based $\ell _1$ regularization family has a significant performance improvement compared to that of their original form which is consistent with the findings in the source localization simulation above, and this phenomenon can be observed in Fig. 5, where the reconstruction results of the GFT-based $\ell _1$ regularization family method are significantly better than the original method, especially the reconstruction of the source signal trend during some time intervals.

3.2 Real data experiments

We further validated the proposed methodology on a real dataset that is publicly accessible through the MNE-Python package [38]. In this dataset acquirement, checkerboard patterns were presented into the left and right visual field, interspersed by tones to the left or right ear with a stimuli interval of 750 ms. The subject was asked to press a key with the right index finger as soon as possible after the appearance of a smiley face was presented at the center of the visual field [40]. Evoked Response Potentials (ERP) were extracted from the MEG measurements, and then we averaged these ERPs for source reconstruction under MNE, MCE, $\ell _{2,1}$, dSPM, sLORETA with and without the proposed GFT operations. The averaged spikes are shown in Fig. 6, and the reconstructed source distributions are shown in Fig. 7.

Fig. 6

Averaged MEG time series plot and topographies

Fig. 7

Reconstructed source activation patterns from MEG data

From Fig. 7, we can see that the source area estimated by MNE, MCE, $\ell _{2,1}$, sLORETA, and dSPM is highly broad. By contrast, the proposed GFT-based methods provide more sparse focal source reconstructions. Moreover, the reconstructed focal for the proposed GFT-based methods falls primarily on areas with the strongest source signal, while others would spread to several regions. Obviously, the spatial graph filter in the proposed GFT-based methods promotes a concentrated and accurate estimation of the visual zone.

4 Discussion

We showed that the subspace composed of low-frequency spatial graph filters can be used to approximate the compact and extended source activation patterns, which is a non-parametric approach, thus empowering the classical ESI methods to have better reconstruction performance for source extents reconstruction. This approach is an easy and effective pre-processing technique compared to predefined spatial gradient based methods or data-driven methods in a deep learning paradigm. While we validated 5 classical methods, it is worth noting that most of the other classical approaches, such as recursive multiple signal classifier (MUSIC), recursively applied and projected MUSIC (RAP MUSIC), low-resolution brain electromagnetic tomography (LORETA), FOCUSS, weighted minimum norm (WMN), shrinking LORETA-FOCUSS, hybrid weighted minimum norm method, recursive sLORETA-FOCUSS, standardized shrinking LORETA-FOCUSS (SSLOFO), etc., can use the proposed GFT in the source space before applying each specific method.

5 Conclusion

In this study, we improved the classical source localization methods by using spatial graph filters to solve the inverse problem of ESI. The proposed methodology enjoys the advantage of reconstructing focal source extents with sparsity and minimizing the impact of the noise by transforming the estimation of the source signal into a projected subspace spanned by spatial frequency graph filters. Numerical experiments demonstrated that the proposed method performs particularly well on source extents, yields excellent robustness when the SNR level is low, and can better reconstruct the source time series, specifically the performance of the $\ell _1$ family regularization can be greatly improved with the GFT. In the experiment on real data we performed, the proposed methodology provides a satisfactory reconstruction with more concentrated source distribution and more stability to noise than benchmark algorithms.