Source localization using recursively applied and projected MUSIC with flexible extent estimation

Magneto- and electroencephalography (M/EEG) are widespread techniques to measure neural activity in-vivo at a high temporal resolution but low spatial resolution. Locating the neural sources underlying the M/EEG poses an inverse problem, which is ill-posed. We developed a new method based on Recursive Application of Multiple Signal Classification (MUSIC). Our proposed method is able to recover not only the locations but, in contrast to other inverse solutions, also the extent of active brain regions flexibly (FLEX-MUSIC). This is achieved by allowing it to search not only for single dipoles but also dipole clusters of increasing extent to find the best fit during each recursion. FLEX-MUSIC achieved the highest accuracy for both single dipole and extended sources compared to all other methods tested. Remarkably, FLEX-MUSIC was capable to accurately estimate the level of sparsity in the source space (r = 0.82), whereas all other approaches tested failed to do so (r ≤ 0.18). The average computation time of FLEX-MUSIC was considerably lower compared to a popular Bayesian approach and comparable to that of another recursive MUSIC approach and eLORETA. FLEX-MUSIC produces only few errors and was capable to reliably estimate the extent of sources. The accuracy and low computation time of FLEX-MUSIC renders it an improved technique to solve M/EEG inverse problems both in neuroscience research and potentially in pre-surgery diagnostic in epilepsy.

1. Introduction

In this paper, we present a novel approach for solving the inverse problem of magneto and electroencephalography (M/EEG) using truncated recursively applied multi-signal classification for sources with variable coherent (FLEX-MUSIC). The EEG and MEG inverse problem is a fundamental challenge in the field of neuroscience, as it involves inferring the underlying neural activity that generates a given set of EEG or MEG measurements (Nunez and Srinivasan, 2006; He et al., 2018; Awan et al., 2019; Michel and Brunet, 2019). The problem is that many different configurations of brain activity can cause the same signal measured on the scalp. Traditional methods for solving this inverse problem rely on mathematical assumptions that are often not aligned with biophysical models of the brain. Others rely on statistical techniques or optimization algorithms, which can be computationally expensive and may not scale well to large datasets. Despite the underdetermined nature of the M/EEG inverse problem, inverse solutions help researchers and clinicians alike. M/EEG inverse solutions help researchers gain insights into the spatio-temporal workings of the brain during, e.g., perceptual and/ or cognitive processing (Feige et al., 2005; Luck, 2014; Kornmeier et al., 2019). M/EEG inverse solutions allow for the identification of functional brain networks and is also used in neurofeedback applications where participants learn to control brain activity in different frequency bands within a specified location (van Lutterveld et al., 2017). Clinicians utilize M/EEG inverse solutions to categorize brain abnormalities or to guide brain surgery, e.g., in drug-resistant epileptic patients (Ebersole, 1994; Lantz et al., 1996; Aydin et al., 2015; Willemse et al., 2016; Sharma et al., 2018).

There are several existing approaches for solving the M/EEG inverse problem (henceforth referred to as solvers), each with its own strengths and limitations. A popular class of methods are the minimum norm estimates (MNE, Hamalainen, 1984; Hämäläinen and Ilmoniemi, 1994; Pascual-Marqui, 1999) that aim to find a solution that explains the observed EEG data with minimal energy of the sources. While these approaches typically incorporate the L2-norm of the source, L1-type solvers have been proposed and termed minimum current estimates (MCE) in the domain of M/EEG inverse problems (Beck and Teboulle, 2009). These L1-type solvers find sources that are sparse in nature since the imposed L1-penalty ultimately sets most dipole values close to zero.

Low-resolution tomography (LORETA) and its iterations, standardized and exact LORETA, fall within this class of solvers (Pascual-Marqui, 1999, 2002, 2007). While these approaches are fast to compute and easy to interpret, these approaches often produce blurred solutions, despite their ability to correctly localize source maxima.

Beamformers are a class of spatial filtering techniques that are commonly used to solve the M/EEG inverse problem. These methods involve constructing a set of spatial filters that are applied to the measured signals to estimate the neural activity at each location in the brain. Beamformers can be effective in certain scenarios but struggle to correctly localize multiple correlated sources or spatially coherent source activity. Furthermore, minimum variance beamformers are known to be sensitive to errors in the forward model. Robust minimum variance beamformers have shown to mitigate the impact of modelling errors recently (Hosseini et al., 2018). A commonly used proponent is the linearly constrained minimum variance (LCMV) Beamformers (Van Veen et al., 1997; Grech et al., 2008). A novelty in the field of Beamformers are the multiple constrained minimum variance (MCMV) Beamformers, which alleviate the problem of correlated sources (Nunes et al., 2020).

Bayesian methods provide a framework for incorporating prior knowledge and uncertainty into the inverse solving process (Friston et al., 2008; Grech et al., 2008; Wipf and Nagarajan, 2009; Wipf et al., 2010). These approaches involve constructing a probabilistic model of the forward and inverse problems, and then apply Bayesian inference to estimate the posterior distribution of the neural activity given the measured M/EEG signals. Bayesian approaches, like sparse Bayesian learning (SBL, Friston et al., 2008; Wipf and Nagarajan, 2009), are capable to accurately localize sources under different sparsity assumptions. However, Bayesian optimization can often be computationally intensive. This problem is amplified when many dipoles are present in the source space, exacerbating the computational load.

A novel class of inverse solvers arose in the past decade that utilize the recent advances in machine learning, predominantly artificial neural networks (ANNs) to solve M/EEG inverse problems. These approaches require training an ANN to produce a source estimate based on simulated pairs of source and M/EEG activity and achieve high accuracy compared to many conventional methods (Cui et al., 2019; Hecker et al., 2020, 2022; Pantazis and Adler, 2021). ANNs are prone to biases in the training data, wherefore their application is yet limited.

Multiple Signal Classification (MUSIC, Mosher and Leahy, 1998) and recursively applied (RAP-) MUSIC (Mosher and Leahy, 1999) are popular approaches for solving the M/EEG inverse problem. Both methods are based on the concept of subspace estimation, which involves estimating the subspace of the neural activity from the measured M/EEG signals using singular value decomposition (SVD). While MUSIC calculates an inverse solution by applying an SVD on the data covariance matrix in a single step, RAP-MUSIC extends the MUSIC method by repeatedly applying the MUSIC algorithm to the measured M/EEG signals. At each iteration, the estimated signal subspace is used to update the estimate of the neural activity, and the updated estimate is then used to update the estimate of the noise subspace. This recursive process continues until convergence is achieved. An improvement of the RAP-MUSIC algorithm was proposed by truncating the recursively calculated subspace with each iteration (TRAP-MUSIC, Mäkelä et al., 2018). This effectively removes residual variance that could not be explained in the prior iterations which leads to disturbances in localization (“RAP dilemma”, cf. Figure 2 in Mäkelä et al., 2018).

Algorithms that follow the RAP-MUSIC-scheme have in common that a spatially coherent source patch will not be detected reliably. Approaches to find multiple dipoles per recursion were proposed by Mosher and Leahy (1998) and Katyal and Schimpf (2004), albeit with factorial increase in computational complexity. Liu and Schimpf (2006) introduced a more computationally efficient way to find extended source clusters by applying computing a weighted Minimum Norm Estimate (wMNE) solution on all estimated single dipoles and their respective neighbors. This step, however, can introduce (1) spurious sources at locations in which a single dipole would have explained sufficient parts of the signal subspace and (2) the approach is limited to finding sources of larger extent beyond single dipoles and their first-order neighbors.

In this paper we summarize the results of our endeavor to overcome this limitation in creating an analysis algorithm that is capable to effectively solve the inverse problem in constellations where established methods do have limitations.

2. Methods

The M/EEG inverse problem refers to the process of inferring the underlying neural activity that generates a given set of EEG measurements. This problem is a fundamental challenge in the field of neuroscience, as it involves understanding the spatio-temporal patterns of neural activity that underlie various cognitive and behavioral processes.

The EEG inverse problem is described as correctly identifying the source matrix J ∈ ℝ^p×t that, when multiplied by the leadfield matrix L ∈ ℝ^q×p (also referred to as gain-matrix) produces the observed EEG data matrix M ∈ ℝ^q×t with q electrodes, p dipoles (i.e., positions in the brain) and t time points.

The propagation of neural currents J through the leadfield L is thus defined as:

$\begin{array}{ll}M=LJ+\zeta ,& (1)\end{array}$

where ζ is the noise in the EEG data. For simplicity, we assume dipoles with fixed orientation perpendicular to the cortical surface.

2.1. Multi-Signal Classification (MUSIC) and its iterations

MUSIC approaches can be interpreted as algorithms that aim to select candidate dipoles in the brain that explain the signal of the EEG data. This is accomplished by calculating the signal and noise subspace of the data covariance C ∈ ℝ^q×q:

$\begin{array}{ll}\begin{array}{c}C=M{M}^T\\ U_s{D}_s{U}_s^T=C_0,\end{array}& (2)\end{array}$

where M^T denote the transpose of the EEG data matrix M and C₀ denotes the covariance matrix of the noiseless data. U_s denotes the eigenvectors and D_s the eigenvalues that both belong to the signal subspace. Signal and noise are sought to be disentangled by selecting only the first n eigenvalues of the covariance matrix UDU^T = C. While this selection is inherently difficult, it was recommended to overestimate n to avoid losing parts of the signal subspace (Mosher and Leahy, 1999). We followed a different approach by algorithmically selecting the set of eigenvalues belonging to the signal subspace. First, eigenvalues were normalized to a maximum of 1 by dividing all eigenvalues by the largest eigenvalue, yielding the normalized set of eigenvalues $\tilde{D}$. We then calculate the difference from each eigenvalue to the next eigenvalue, yielding ${\tilde{\text{d}}}_{\text{Δ}}$. Let ${\tilde{\text{d}}}_{\text{Δ}}$ be the set of eigenvalues and ϵ = 0.01 be the relative selection criterion, which was determined empirically during our testing phase. The smallest eigenvector ñ that belongs to the signal subspace is defined as

$\begin{array}{ll}\begin{array}{l}\begin{array}{c}UD{V}^T=C\\ d=[D_{1,1},D_{2,2},\mathrm{..}{D}_{q,q}]\\ {\tilde{d}}_i=\frac{d_i}{max(d)}\\ {\tilde{d}}_{\text{Δ},i}={\tilde{d}}_{i+1}-{\tilde{d}}_i,\forall i\in \{2,\dots ,q\}\\ ñ=i|{\tilde{d}}_i<ϵ\le {\tilde{d}}_{i-1}\end{array}\hfill \end{array}& (3)\end{array}$

The estimated signal subspace is thus U_s = U(1:ñ) and the projection to the signal space is defined as $P_s=U_s{U}_s^T$.

The MUSIC localizer is then calculated as follows:

$\begin{array}{ll}\begin{array}{l}{\mu }_p=\frac{\left|\right|P_{s}l(p)|{|}^2}{\left|\right|l(p)|{|}^2},\hfill \end{array}& (4)\end{array}$

where l(p) denotes the p^th column of the leadfield matrix L. The resulting localizer μ is finally filtered to only contain values above a certain criterion (typically between 0.9 to 0.99).

Recursively applied MUSIC (RAP-MUSIC) also makes use of the signal subspace and its projection by iteratively selecting candidate dipoles as follows. We henceforth describe the RAP-MUSIC algorithm.

The first candidate is the dipole with the largest source amplitude in the MUSIC localizer as described above (Eq. 4). Let Î₁ be the topography of the initially selected candidate at iteration i we construct a set of topographies B = [Î₁, ..., Î_i] that stores all topographies of the selected candidates. Using the set of topographies B ∈ ℝ^q×i at iteration i we define the out-projector matrix Q_i:

$\begin{array}{ll}{Q}_i=I-B_i{B}_i^{†},& (5)\end{array}$

where I ∈ ℝ^q×q denotes the identity matrix and B^† denotes the Moore-Penrose pseudo inverse of B.

The updated covariance matrix C_i is then calculated by multiplying the out-projector matrix by the signal subspace and the new signal subspace is thusly calculated:

$\begin{array}{ll}\begin{array}{c}C=Q{U}_i\\ C=U_i{D}_i{U}_i^T\end{array}& (6)\end{array}$

We then calculate the new signal subspace projection as $P_i=U_i(1:ñ)U_i{(1:ñ)}^T$. Note, that the number of selected components can be truncated to $U_i(1:n-\tilde{k}-1)$ to alleviate the problem of the RAP-dilemma (Mäkelä et al., 2018). The (T)RAP-MUSIC localizer is then calculated as follows:

$\begin{array}{ll}\begin{array}{l}{\mu }_i(p)=\frac{\left|\right|P_i{Q}_{i}l(p)|{|}^2}{\left|\right|Q_{i}l(p)|{|}^2},\hfill \end{array}& (7)\end{array}$

where l(p) denotes the p^th column of the leadfield matrix L. The dipole at which u_i is maximal is selected as the new candidate c_i. The equations 5, 6 & 7 are iterated while i < q or a stopping criterion is met. The stopping criterion typically is met when the maximum of the (T)RAP-MUSIC localizer u_i falls below a threshold (in our case 0.975), i.e., when no dipole is able to explain the signal subspace projection sufficiently.

2.2. FLEX-MUSIC

A new dipole is selected in each iteration of the RAP-MUSIC algorithm based on the current subspace projection P. As described by Mäkelä et al. (2018), the selection of the most optimal dipoles often leaves some residual of the signal subspace projection to be explained. One reason for this is that neural sources can often not be represented by a single dipole due to the functional coherence in the cortex that is reflected by locally smooth activity of varying extent. To overcome this limitation we extended the set of candidate dipoles by multiple sets of smoothly distributed dipole clusters.

We have therefore calculated gradients $G_k\in {ℝ}^{p\times p}$ for increasing smoothness orders k ∈ 1, 2, ..., 8. Each gradient G_k transforms the original leadfield matrix L to leadfield matrices of increasing smoothness orders L_k as based on the adjacency matrix A ∈ ℝ^q×q:

$\begin{array}{ll}\begin{array}{l}\begin{array}{c}{G}_1=I\\ G_2=laplacian(A)\\ G_k=G_{k-1}A|k>2\\ L_k=L{G}_k,\\ \end{array}\hfill \end{array}& (8)\end{array}$

where S is the largest order of smoothness and I ∈ ℝ^p×p is the identity matrix. In summary, Eq. 8 states that we create a set leadfield matrices G_k with progressing smoothness. This is achieved using the Laplacian of the adjacency matrix. The adjacency matrix, also referred to as neighborhood matrix describes which dipoles in our model are directly connected (i.e., are neighbors).

The selection of the highest smoothness order k depends on the upper boundary of smoothness to be assumed in the source model and should be adjusted depending on the number of dipoles in the source model. Note, that the original leadfield, denoted as L₁, remains, since the respective gradient G₁ is the identity matrix.

We calculate the FLEX-MUSIC localizer at iteration i by

$\begin{array}{ll}\begin{array}{l}{\mu }_{i,k}=\frac{\left|\right|P_i{Q}_i{L}_k|{|}^2}{\left|\right|Q_i{L}_k|{|}^2},\hfill \end{array}& (9)\end{array}$

yielding one localizer for each smoothness order k.

The dipole or dipole cluster $\widehat{p}$ (depending on the estimated optimal smoothness order $\widehat{k}$) at which μ_{i, k} is maximal is selected as the new candidate.

We then update the set of topographies B with the newly added topography $L_{\widehat{k},\widehat{p}}$. Furthermore, we update the source covariance matrix S by adding the column vector of $G_{\widehat{k},\widehat{p}}$:

$\begin{array}{ll}{S}_i=S_{i-1}+G_{\widehat{k},\widehat{p}}& (10)\end{array}$

FLEX-MUSIC iterates the equations 5, 6, 8, 9 and 10. No truncation, as done in TRAP-MUSIC, was applied since that yielded better results during testing.

Unlike MUSIC, the recursive approaches (e.g., RAP-, TRAP- and FLEX-MUSIC) require a final estimation of the current source density Ĵ ∈ ℝ^p×t after candidate selection. We use the source covariance matrix S to calculate a weighted minimum-norm-like solution:

$\begin{array}{ll}Ĵ=S{L}^T{(LS{L}^T)}^{-1}M,& (11)\end{array}$

where M ∈ ℝ^q×t is the EEG or MEG data matrix.

In summary, FLEX-MUSIC further alleviates the RAP-dilemma of residual variance in the signal subspace projector P. This is achieved by adding clusters of neighboring dipoles to the set of candidate dipoles at each recursion step. The algorithm is still capable to localize single dipoles since they remain part of the set of candidates (cf. Eq. 8). This renders FLEX-MUSIC a flexible solution to the M/EEG inverse problem in which the spatial extent of neural activations is often unknown a priori. In this way, FLEX-MUSIC increases the probability to identify extended sources and integrate their extensions into the result source space, instead of either ignoring the extension or treating the neighboring dipoles as separate sources.

2.3. Evaluation

In order to evaluate the proposed method, we simulated pairs of source- and EEG- data using an anatomical template brain “fsaverage” (Fischl et al., 1999) by the Freesurfer image analysis suite¹. EEG simulations were carried out using a precomputed three shell boundary element method (BEM; Fuchs et al., 2002) forward solution as provided by mne-python (v20.3, Gramfort et al., 2013). Each shell (brain, skull & scalp tissue) was composed of 5120 vertices. The conductivity was set to 0.3S/m² for brain and scalp tissue, and 0.06S/m² for the skull.

The source model was chosen with p = 1, 284 dipoles with icosahedral spacing. For the EEG electrodes we used the Biosemi 64-channel layout consisting of q = 64 electrodes of the 10-20 system. Using the forward model and the parameters described, we calculated a leadfield L ∈ ℝ^q×p.

We evaluate our proposed method FLEX-MUSIC by comparing it to a diverse set of other inverse algorithms including TRAP-MUSIC, eLORETA, a sparse Bayesian learning (SBL) approach called Convexity Champagne and the Multiple Constrained Minimum Variance Beamformer (MCMV, Mosher and Leahy, 1999; Pascual-Marqui, 2007; Wipf and Nagarajan, 2009; Mäkelä et al., 2018; Nunes et al., 2020; Cai et al., 2022).

Motivation is given for the choice of each method for solving the EEG inverse problem. TRAP-MUSIC was chosen as one of the latest developments of the recursive MUSIC approaches. eLORETA is a popular choice in many EEG studies with theoretically low localization errors, rendering it the most suitable choice within the minimum-norm family. Convexity Champagne was chosen as a very recent improvement to the Champagne algorithm. Champagne was shown to produce fast and accurate solutions within the framework of empirical Bayes (Wipf et al., 2010). MCMV is a similar approach as the linearly constrained minimum variance (LCMV) beamformer. However, it is designed to be less prone to correlations between sources, rendering it a useful innovation over LCMV. All methods were implemented in-house and are available in our python package invertmeeg².

Optimal regularization of eLORETA, MCMV and Convexity Champagne was achieved using generalized cross validation (GCV, Grech et al., 2008) on a set of 7 regularization parameters λ = {10⁻³..., 10³}. The MUSIC-type methods (FLEX- and TRAP-MUSIC) do not require the regularization parameters since the noise is estimated by selecting the signal subspace.

A set of 1,000 samples consisting of ground truth sources J ∈ ℝ^p×t and corresponding EEG M ∈ ℝ^p×t was simulated in order to evaluate the accuracy of all solvers.

The number of simulated consecutive time points was set to t = 20. The simulation parameters are outlined in Table 1. Half of all samples contained single dipoles, whereby the other half contained samples of extended dipole clusters with coherent activity over time. The cluster size was varied in terms of neighborhood orders, whereas an order of 1 indicates a single dipole and an order of 2 indicates a dipole including all its neighbors. The source time course was generated as random sequence of a colored frequency spectrum as described by the $P(f)=\frac{1}{f^{beta}}$, where beta controls the level of temporal smoothness. Noise was generated as random spatio-temporal white noise with inter-channel correlation between -1 and 1. The noise was added to the EEG matrix such that a random signal-to-noise ratio (SNR) within the given range outlined in Table 1 was achieved. Only white noise is considered since the presence of colored noise in real data can be handled by whitening the EEG data as a preprocessing step.

TABLE 1

Table 1. Simulation parameters.

We calculated inverse solutions to each of the simulated samples of EEG data using the different methods, as described above. Accuracy of the individual inverse solutions is quantified by calculating the mean localization error (MLE), Earth Mover's Distance (EMD, Hitchcock, 1941) and the mean squared error (MSE). Furthermore, we quantified the sparsity of each inverse solution.

We calculated the MLE by first identifying the local maxima of the ground truth source matrix J and the inverse solution Ĵ and then calculating the minimum Euclidean distance between each true dipole location and all estimated dipole locations.

MSE quantifies how close the estimated dipole moments (in nAm) are to the true dipole moments and was calculated as follows:

$\begin{array}{ll}\begin{array}{l}MSE(J,\widehat{J})=\left|\right|J-\widehat{J}|{|}_2\hfill \end{array}& (12)\end{array}$

EMD is a measure that calculates the distance between two distributions. It is a suitable method to quantify the accuracy of an inverse solution since, unlike MSE, it takes into account the distance between dipole locations. It was calculated as follows:

$\begin{array}{ll}EMD(J,Ĵ)={\displaystyle \sum _p}D(J-Ĵ),& (13)\end{array}$

where D ∈ ℝ^p×p is the distance matrix containing the Euclidean distance between each dipole pair. Prior to calculating the EMD, we have computed the absolute mean over time points for J and Ĵ and normalized them.

Sparsity was calculated by first normalizing the columns of the estimated source matrix Ĵ ∈ ℝ^p×t to unit length by division of the respective columns L2-norm. The L1 norm of the normalized matrix ||Ĵ||₁ was then calculated, yielding a metric with an inverse relationship to sparsity. According to the dogma of Occam's Razor, which states that complexity should not be posited without necessity, we can assume that sparse solutions make the fewest assumptions about the brain's activity and are therefore preferred, given that they explain a sufficient amount of the data.

In summary, the evaluation metrics described above capture the accuracy of estimated local maxima positions (MLE), the accuracy of the global pattern of the inverse solution (EMD), the accuracy of dipole moments (MSE) and the sparsity of the solution (L1 norm).

3. Results

We calculated the accuracy of all inverse algorithms as described in the previous section. Exemplary samples of ground truth source activity and estimated sources are shown in Figure 1.

FIGURE 1

Figure 1. Examples of ground truth and estimated sources. Exemplary plots of ground truth sources and the estimated sources of FLEX-MUSIC and comparative approaches. For demonstration purposes, samples were selected based on visibility on the right lateral view. (First Row) Single non-extended source. (Seconds Row) Multiple non-extended sources. (Third Row) Single extended source. (Fourth row) Multiple extended sources. Colorbars were adjusted for eLORETA and MCMV to improve visibility of the source pattern. FLEX-MUSIC visibly recovers the actual source extent whereas all other approaches tested show biases toward single dipoles or extended dipoles.

Figure 2 depicts the Mean Localization Error (MLE), Earth Mover's Distance (EMD) and the Mean Squared Error (MSE) for all inverse solutions of each solver. All metrics are summarized in Table 2. Samples were divided into those containing single dipole sources and those containing extended source clusters. FLEX- and TRAP-MUSIC show overall lowest MLE and EMD for single dipole sources when compared to all other solvers. Notably, the median MLE is zero for both MUSIC-based methods, and there was no significant difference in MLE for single-dipole sources (p = 0.96, t = 0.36, d = 0.02). FLEX- and TRAP-MUSIC did also not differ significantly in EMD (p = 0.95, t = 0.73, d = 0.05) and MSE (p = 1.00, t = 0.31, d = 0.02). Champagne produced the next best accuracies, eLORETA and the MCMV Beamformer show comparatively poor accuracy in correctly estimating the source distribution as depicted in relatively high EMD and MSE, whereas the Convexity Champagne solver lies in-between.

FIGURE 2

Figure 2. Evaluation of all solvers. Boxplots depict the accuracy of FLEX-MUSIC and all other solvers tested. (Left) Mean Localization Error in mm. (Center) Earth Mover's Distance (EMD), (Right) Mean Squared Error. Blue: Single dipoles. Orange: Extended dipoles. Note, that FLEX-MUSIC achieves competitive accuracy for single-dipole sources and the highest accuracy for samples containing spatially extended sources.

TABLE 2

Table 2. Medians of all metrics.

The advantage of the proposed FLEX-MUSIC solver becomes most apparent for extended sources. While TRAP-MUSIC fails to accurately localize sources with spatial extent, FLEX-MUSIC retains the lowest MLE and EMD compared to all other solvers (Figure 2).

Figure 3 shows the sparsity of the produced inverse solutions of all solvers, defined as the L1-norm of the L2-normalized source estimate Ĵ. We find that FLEX- and TRAP-MUSIC exhibit the highest sparsity, MCMV and eLORETA the lowest, and Convexity Champagne lies in the middle. Interestingly, only FLEX-MUSIC shows a clear difference in sparsity between samples containing single dipoles and those containing extended dipole clusters, which is consequence of its flexibility to estimate sources of varying extent.

FIGURE 3

Figure 3. Sparsity of inverse solutions. Boxplot depicts the L1 norm for all inverse solvers tested. The L1 norm reflects the level of “non-sparsity”.

This aspect is shown in detail in Figure 4, depicting each solver's capability to recover the spatial level of sparsity in the ground truth. FLEX-MUSIC shows the highest correlation (r = 0.82) between the sparsity in the ground truth sources and the sparsity in the predicted sources. TRAP-MUSIC shows a high correlation for highly sparse samples and is biased for less sparse samples. Convexity Champagne produced solutions that were often less sparse that the ground truth, whereas MCMV and eLORETA had a strong bias toward finding less sparse activations. Despite the strong bias, the sparsity of eLORETA was significantly correlated with the sparsity in the ground truth samples.

FIGURE 4

Figure 4. Accuracy in extent estimation. Scatter plots depict the relationship between the L1 norm in the ground truth and the L1 norm in the prediction. The L1 norm reflects the level of “non-sparsity.” Note, that only FLEX-MUSIC is capable to reproduce the true sparsity with high correlation (r = 0.82). ***p < 0.001.

Next, we tested the dependence of the inverse solution accuracy on varying levels of noise. For this comparison we have combined samples containing single dipoles and those containing extended dipole clusters (Figure 5). The graph shows that FLEX-MUSIC yields inverses solutions with the highest accuracy regardless of the level of noise in the EEG data.

FIGURE 5

Figure 5. Dependence on SNR. Accuracy of all solvers separated by the signal-to-noise ratio (SNR) in the simulated EEG samples. (Left) Mean localization Error for different levels of noise in the EEG data. (Center) Earth Mover's Distance for different levels of noise in the EEG data. (Right) Mean Squared Error for different levels of noise in the EEG data. Note that the advantage of FLEX-MUSIC over other solvers persists for varying levels of noise.

Finally, we analyzed the accuracy of all solvers depending on a varying number of dipoles/ dipole clusters in the ground truth (Figure 6). FLEX-MUSIC achieves again the highest accuracy on all metrics regardless of how many dipoles/ dipole clusters were present in the ground truth source. Similar to Champagne, FLEX-MUSIC showed remarkably low MSE when multiple sources were active simultaneously.

FIGURE 6

Figure 6. Dependence on number of active sources. Accuracy of all solvers for varying numbers of active sources within the simulated EEG samples. Error bars depict standard errors of the mean (SEM).

3.1. Computational expense

The computational expense to calculate the EEG inverse operators is presented in Figure 7 (right). The framework for these calculations was the invertmeeg library which was developed in-house in python (https://github.com/lukethehecker/invert). Although FLEX-MUSIC has some additional processing steps compared to TRAP-MUSIC, we find that the median computation time differs only slightly (Δt = 0.042s, p = 1.72·10⁻⁶, t = 5.25, d = 0.23). eLORETA required similar computation times of 0.32s. Convexity Champagne required the longest computation time of 2.66s. Note, that eLORETA, Convexity Champagne and MCMV were re-computed 7 times for varying levels of regularization. If the optimal regularization parameter was known in advance, the computation time could be reduced seven-fold.

FIGURE 7

Figure 7. Time of computation. Boxplot depicts the computation time in seconds for all inverse solvers tested. Convexity Champagne requires ten times longer computation time compared to FLEX-MUSIC.

4. Discussion

In this work we have presented FLEX-MUSIC, a new approach to solve the M/EEG inverse problem embedded in the recursive MUSIC scheme. Similar to the well-known RAP-MUSIC techniques, it iteratively adds candidate positions to the set of active sources. In addition to the RAP-MUSIC approaches it adds an extended dictionary of increasingly smooth source patches.

We have shown that FLEX-MUSIC works as well as TRAP-MUSIC in scenarios where the EEG was produced by singular dipole sources. Furthermore, we have shown that FLEX-MUSIC accurately estimates the spatial extent of the underlying sources. Strikingly, of all inverse solvers tested, only FLEX-MUSIC was capable of estimating the level of sparsity in the ground truth, showing the highest correlation between the sparsity levels in the ground truth sources and the estimated sources or r = 0.82, whereas the next best solver achieved only a correlation of r = 0.18. We consider this aspect to be the most compelling argument for FLEX-MUSIC, since the correct estimation of the extent of neural generators underlying an M/EEG signal is of high interest for multiple reasons.

First, as stated above, RAP-MUSIC approaches suffer in general from what is called the RAP dilemma, i.e., the interference of residual variance from previous iterations. Since FLEX-MUSIC is capable to explain a larger portion of the subspace, it effectively diminishes the residual variance from the previous iterations, leading to more accurate inverse solutions throughout all iterations.

Second, good estimations of the extent of neural sources is of high interest in the pre-surgical diagnostic in epilepsy. Not only does FLEX-MUSIC more reliably find the true location of the source maxima regardless of their extent (cf. Figure 2), it is also capable to isolate the prospective resection area in the brain.

As stated earlier, the recursive MUSIC approaches are a suitable option to solve inverse problems where the solution space is considerably large due to their low computation time and built-in regularization. While the presented Convexity Champagne has shown competitive accuracy in many cases, we showed that the computational expense was almost 10 times higher in our setting. We further argue that the difference in computational expense may further increase with larger source models (i.e., higher number of dipoles), rendering the computation of Bayesian inverse solutions unfeasible fir certain clinical applications.

The simulations presented in this work contained spherical sources, similar as to what FLEX-MUSIC is able to model by progressively smoothing the leadfield. This may have led to optimistic results, and it is expected that deviant source shapes (e.g., elliptical) could cause FLEX-MUSIC to perform worse. However, we have made efforts to simulate a broad set of neural activity with up to 10 simultaneously active sources. Diverse spatial patterns of neural activity form when many spherical sources are combined, e.g., through overlap. FLEX-MUSIC performed well even under these conditions, as was shown in Figure 6, which indicates robustness of the algorithm.

A potential weakness of our proposed FLEX-MUSIC algorithm is that mesoscale brain activity may not be sufficiently modelled with single dipoles and smooth dipole clusters. Various shapes, e.g., elliptical coherent sources, may be involved in real-world M/EEG recordings of brain activity. However, we expect that for a sufficiently large source model, FLEX-MUSIC should still be able to reconstruct deviant shapes of sources using the circular smooth patches. Future improvements of the FLEX-MUSIC algorithm could involve changing the way the dictionary of candidate dipoles or dipole clusters is applied. One potentially fruitful modification may be to find the set of active dipoles that explain the signal subspace in each iteration is found in a data-driven manner, unlike the dictionary-like search that we proposed in the present work. Another idea would be a combination of SBL and FLEX-MUSIC, in which a small set of active dipoles are added to the source covariance matrix (cf. Eq. 10) based on the maximum a posteriori (MAP) estimation. An interesting development of the recursive MUSIC algorithms was presented recently by Adler et al. (2022), called Alternating Projections (AP). The approach yielded lower MLE compared to RAP- and TRAP-MUSIC, especially under low SNR conditions. It may be promising to translate the idea of FLEX-MUSIC to the domain of AP, which could potentially further increase the accuracy of the AP-approach in scenarios where spatially coherent sources are expected.

A very important finding of the present work, beyond the performance benefits of FLEX-MUSIC, is that different inverse solutions applied to the same data set can produce quite divergent results. In the case of real EEG data, where a ground truth is not available, we strongly recommend to use different methods for EEG source calculations in parallel and compare the results with each other. As we have shown in this work, combining FLEX-MUSIC and a Champagne algorithm may be useful, although Champagne may overestimate the size of neural sources in some cases.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/LukeTheHecker/flex_paper_analyses.

Author contributions

LH is the author of the idea and implemented the inverse solution and wrote the initial manuscript. LH, JK, and LT designed the evaluation procedure. All authors reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The project was funded by the European Campus (Eucor) seeding money.

Acknowledgments

We acknowledge support by the Open Access Publication Fund of the University of Freiburg.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

1. ^http://surfer.nmr.mgh.harvard.edu/

2. ^https://github.com/lukethehecker/invert

References

Adler, A., Wax, M., and Pantazis, D. (2022). “Brain source localization by alternating projection,” in 2022 IEEE 19th International Symposium on Biomedical Imaging (ISBI) (IEEE), 1–5. doi: 10.1109/ISBI52829.2022.9761604

CrossRef Full Text | Google Scholar

Awan, F. G., Saleem, O., and Kiran, A. (2019). Recent trends and advances in solving the inverse problem for EEG source localization. Inverse Probl. Sci. Eng. 27, 1521–1536. doi: 10.1080/17415977.2018.1490279

CrossRef Full Text | Google Scholar

Aydin, Ü., Vorwerk, J., Dümpelmann, M., Küpper, P., Kugel, H., Heers, M., et al. (2015). Combined EEG/MEG can outperform single modality EEG or MEG source reconstruction in presurgical epilepsy diagnosis. PLoS ONE 10, e0118753. doi: 10.1371/journal.pone.0118753

PubMed Abstract | CrossRef Full Text | Google Scholar

Beck, A., and Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imaging Sci. 2, 183–202. doi: 10.1137/080716542

CrossRef Full Text | Google Scholar

Cai, C., Kang, H., Hashemi, A., Chen, D., Diwakar, M., Haufe, S., et al. (2022). Bayesian algorithms for joint estimation of brain activity and noise in electromagnetic imaging. IEEE Trans. Med. Imaging. 42, 762–773. doi: 10.1109/TMI.2022.3218074

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, S., Duan, L., Gong, B., Qiao, Y., Xu, F., Chen, J., et al. (2019). EEG source localization using spatio-temporal neural network. China Commun. 16, 131–143. doi: 10.23919/JCC.2019.07.011

CrossRef Full Text | Google Scholar

Ebersole, J. S. (1994). Non-invasive localization of the epileptogenic focus by EEG dipole modeling. Acta Neurol. Scand. 89, 20–28. doi: 10.1111/j.1600-0404.1994.tb05179.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Feige, B., Scheffler, K., Esposito, F., Di Salle, F., Hennig, J., and Seifritz, E. (2005). Cortical and subcortical correlates of electroencephalographic alpha rhythm modulation. J. Neurophysiol. 93, 2864–2872. doi: 10.1152/jn.00721.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Fischl, B., Sereno, M. I., Tootell, R. B., and Dale, A. M. (1999). High-resolution intersubject averaging and a coordinate system for the cortical surface. Hum. Brain Mapp. 8, 272–284. doi: 10.1002/(SICI)1097-0193(1999)8:4&lt;272::AID-HBM10&gt;3.0.CO;2-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K., Harrison, L., Daunizeau, J., Kiebel, S., Phillips, C., Trujillo-Barreto, N., et al. (2008). Multiple sparse priors for the M/EEG inverse problem. Neuroimage 39, 1104–1120. doi: 10.1016/j.neuroimage.2007.09.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Fuchs, M., Kastner, J., Wagner, M., Hawes, S., and Ebersole, J. S. (2002). A standardized boundary element method volume conductor model. Clin. Neurophysiol. 113, 702–712. doi: 10.1016/S1388-2457(02)00030-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Gramfort, A., Luessi, M., Larson, E., Engemann, D. A., Strohmeier, D., Brodbeck, C., et al. (2013). MEG and EEG data analysis with MNE-Python. Front. Neurosci. 7, 267. doi: 10.3389/fnins.2013.00267

PubMed Abstract | CrossRef Full Text | Google Scholar

Grech, R., Cassar, T., Muscat, J., Camilleri, K. P., Fabri, S. G., Zervakis, M., et al. (2008). Review on solving the inverse problem in EEG source analysis. J. Neuroeng. Rehabil. 5, 25. doi: 10.1186/1743-0003-5-25

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamalainen, M. (1984). Interpreting Measured Magnetic Fields of the Brain: Estimates of Current Distributions. Univ. Helsinki, Finland Tech. Rep. TKK-F-A559.

Google Scholar

Hämäläinen, M. S., and Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35–42. doi: 10.1007/BF02512476

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B., Sohrabpour, A., Brown, E., and Liu, Z. (2018). Electrophysiological source imaging: a noninvasive window to brain dynamics. Ann. Rev. Biomed. Eng. 20, 171–196. doi: 10.1146/annurev-bioeng-062117-120853

PubMed Abstract | CrossRef Full Text | Google Scholar

Hecker, L., Rupprecht, R., van Elst, L. T., and Kornmeier, J. (2020). ConvDip: a convolutional neural network for better M/EEG Source Imaging. bioRxiv. doi: 10.1101/2020.04.09.033506

CrossRef Full Text | Google Scholar

Hecker, L., Rupprecht, R., van Elst, L. T., and Kornmeier, J. (2022). Long-short term memory networks for electric source imaging with distributed dipole models. Neuroscience. bioRxiv. doi: 10.1101/2022.04.13.488148

CrossRef Full Text | Google Scholar

Hitchcock, F. L. (1941). The distribution of a product from several sources to numerous localities. J. Math. Phys. 20, 224–230. doi: 10.1002/sapm1941201224

PubMed Abstract | CrossRef Full Text | Google Scholar

Hosseini, S. A. H., Sohrabpour, A., Akçakaya, M., and He, B. (2018). Electromagnetic brain source imaging by means of a robust minimum variance beamformer. IEEE Trans. Biomed. Eng. 65, 2365–2374. doi: 10.1109/TBME.2018.2859204

PubMed Abstract | CrossRef Full Text | Google Scholar

Katyal, B., and Schimpf, P. H. (2004). “Multiple current dipole estimation in a realistic head model using R-MUSIC,” in The 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 1 (IEEE), 829–832. doi: 10.1109/IEMBS.2004.1403286

PubMed Abstract | CrossRef Full Text | Google Scholar

Kornmeier, J., Friedel, E., Hecker, L., Schmidt, S., and Wittmann, M. (2019). What happens in the brain of meditators when perception changes but not the stimulus? PLoS ONE 14, e0223843. doi: 10.1371/journal.pone.0223843

PubMed Abstract | CrossRef Full Text | Google Scholar

Lantz, G., Holub, M., Ryding, E., and Rosén, I. (1996). Simultaneous intracranial and extracranial recording of interictal epileptiform activity in patients with drug resistant partial epilepsy: patterns of conduction and results from dipole reconstructions. Electroencephalogr. Clin. Neurophysiol. 99, 69–78. doi: 10.1016/0921-884X(96)95686-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., and Schimpf, P. H. (2006). Efficient localization of synchronous EEG source activities using a modified RAP-MUSIC algorithm. IEEE Trans. Biomed. Eng. 53, 652–661. doi: 10.1109/TBME.2006.870236

PubMed Abstract | CrossRef Full Text | Google Scholar

Luck, S. J. (2014). An Introduction to the Event-Related Potential Technique. MIT Press.

Google Scholar

Mäkelä, N., Stenroos, M., Sarvas, J., and Ilmoniemi, R. J. (2018). Truncated rap-music (trap-music) for MEG and EEG source localization. Neuroimage 167, 73–83. doi: 10.1016/j.neuroimage.2017.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Michel, C. M., and Brunet, D. (2019). EEG Source imaging: a practical review of the analysis steps. Front. Neurol. 10, 325. doi: 10.3389/fneur.2019.00325

PubMed Abstract | CrossRef Full Text | Google Scholar

Mosher, J. C., and Leahy, R. M. (1998). Recursive MUSIC: a framework for EEG and MEG source localization. IEEE Trans. Biomed. Eng. 45, 1342–1354. doi: 10.1109/10.725331

PubMed Abstract | CrossRef Full Text | Google Scholar

Mosher, J. C., and Leahy, R. M. (1999). Source localization using recursively applied and projected (RAP) MUSIC. IEEE Trans. Signal Process. 47, 332–340. doi: 10.1109/78.740118

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunes, A. S., Moiseev, A., Kozhemiako, N., Cheung, T., Ribary, U., and Doesburg, S. M. (2020). Multiple constrained minimum variance beamformer (MCMV) performance in connectivity analyses. Neuroimage 208, 116386. doi: 10.1016/j.neuroimage.2019.116386

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez, P. L., and Srinivasan, R. (2006). Electric Fields of the Brain: The Neurophysics of EEG. Oxford University Press, USA. doi: 10.1093/acprof:oso/9780195050387.001.0001

PubMed Abstract | CrossRef Full Text | Google Scholar

Pantazis, D., and Adler, A. (2021). MEG source localization via deep learning. Sensors 21, 4278. doi: 10.3390/s21134278

PubMed Abstract | CrossRef Full Text | Google Scholar

Pascual-Marqui, R. D. (1999). Review of methods for solving the EEG inverse problem. Int. J. Bioelectromagn. 1, 75–86.

Google Scholar

Pascual-Marqui, R. D. (2002). Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find Exp. Clin. Pharmacol. 24, 5–12.

PubMed Abstract | Google Scholar

Pascual-Marqui, R. D. (2007). Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization. arXiv preprint arXiv:0710.3341.

Google Scholar

Sharma, P., Scherg, M., Pinborg, L. H., Fabricius, M., Rubboli, G., Pedersen, B., et al. (2018). Ictal and interictal electric source imaging in pre-surgical evaluation: a prospective study. Eur. J. Neurol. 25, 1154–1160. doi: 10.1111/ene.13676

PubMed Abstract | CrossRef Full Text | Google Scholar

van Lutterveld, R., Houlihan, S. D., Pal, P., Sacchet, M. D., McFarlane-Blake, C., Patel, P. R., et al. (2017). Source-space EEG neurofeedback links subjective experience with brain activity during effortless awareness meditation. Neuroimage 151, 117–127. doi: 10.1016/j.neuroimage.2016.02.047

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Veen, B. D., van Drongelen, W., Yuchtman, M., and Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44, 867–880. doi: 10.1109/10.623056

PubMed Abstract | CrossRef Full Text | Google Scholar

Willemse, R. B., Hillebrand, A., Ronner, H. E., Peter Vandertop, W., and Stam, C. J. (2016). Magnetoencephalographic study of hand and foot sensorimotor organization in 325 consecutive patients evaluated for tumor or epilepsy surgery. Neuroimage Clin. 10, 46–53. doi: 10.1016/j.nicl.2015.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wipf, D., and Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. Neuroimage 44, 947–966. doi: 10.1016/j.neuroimage.2008.02.059

PubMed Abstract | CrossRef Full Text | Google Scholar

Wipf, D. P., Owen, J. P., Attias, H. T., Sekihara, K., and Nagarajan, S. S. (2010). Robust Bayesian estimation of the location, orientation, and time course of multiple correlated neural sources using MEG. Neuroimage 49, 641–655. doi: 10.1016/j.neuroimage.2009.06.083

PubMed Abstract | CrossRef Full Text | Google Scholar