Recent trends and advances in solving the inverse problem for EEG source localization

ABSTRACT

This paper addresses the recent advancements and trends in the field of electroencephalography (EEG) using inverse problem solutions. Using the EEG data of the brain to gather the information regarding the neuronal current source distribution has been a persisting challenge. Since the EEG inverse problem is ill-posed in nature; therefore, it does not offer a unique result. A trivial and precise solution yields a detailed insight regarding the electrical activity as well as the damaged tissue in the brain. Ordinarily, this problem is solved using the regularization techniques, such as minimum norm estimates, mixed-norm estimate, low-resolution electrical tomography, artificial neural networks, and their modified variants. In this paper, the latest algorithmic developments in solving the EEG inverse problem are reviewed. The optimization rendered by these techniques in accurately solving the neural source localization problem is also discussed. The comparative performance analysis of the recent techniques has been presented. Furthermore, a number of future enhancements have also been proposed to further improve the performance of these state-of-the-art techniques.

KEYWORDS:

• Electroencephalography
• inverse problem
• source localization
• minimum norm estimates
• low-resolution electrical tomography

MSC SUBJECT CLASSIFICATIONS:

• 47A52
• 65F22
• 92C55

1. Introduction

The electroencephalography (EEG) forward problem is defined as the transduction of the neuronal current sources into the scalp electrode potentials. The EEG inverse problem, also referred to as the neural source imaging (NSI) problem, computes the spatial location (or position) of ionic current sources by taking into account the electrode potentials as its inputs [1–4]. The solution to NSI problem provides a detailed spatial–temporal analysis of the brain activity for different functions; therefore, it is fundamental in neuro-engineering. An accurate solution of NSI problem helps in comprehending the inner brain activity, its cognitive processes, and its pathological functioning. Moreover, it helps in identifying the damaged tissues in the brain by estimating conductivity differences [5–7].

Extensive research has been done to accurately solve the inverse problem. The NSI problem does not offer a trivial output due to its ill-posed nature. Hence, in order to find an optimal (or even sub-optimal) solution of this ill-posed problem, the regularization techniques such as low-resolution electromagnetic tomography activity (LORETA), minimum norm estimation (MNE), and their modified variants; namely, weighted-MNE, standardized-LORETA, variable-resolution-electrical-tomography (VARETA), have been widely used [8–10]. Other classical techniques rely upon the utilization of sparsity-inducing penalties to derive an optimal solution. These techniques include minimum current estimate and mixed-norm method (MxNE) [11,12]. Maximum entropy on the mean is another widely used probabilistic technique to solve the inverse problem [13]. These approaches assume a quasi-static approximation in order to mathematically relate the neuronal current sources with the electrode potentials. This mathematical mapping is derived via the lead-field matrix [14]. Techniques involving the brain’s geometric model and application of simplex method to solve the ill-posed problem have been proposed in [15–18]. Other mentionable schemes have been presented in references [19–21].

The aforementioned conventional methods involve convex optimization using L₁ norm [22]. The iterative mixed-norm estimation (irMxNE) combines the iterative reweighting strategy with mixed-norm and quasi-norm. It imposes appropriate penalties on the regions that tend to have lower (or negligible) concentrations of active source [23]. Apart from accurately localizing the sources, this approach also enables the algorithm to compute a precise estimate regarding their strength, in functionally connected regions of the brain.

The electromagnetic properties of the neuronal current sources have been rigorously studied to solve the NSI problem. An effective solving technique via the three-dimensional differential equations has been presented in Afasari et al. [24]. Unlike the previously used integration based techniques, the differential inverse scattering methodology uses eight antennas and differential equations to detect the electric fields rendered by the sources. Finally, Born’s iteration technique and finite element method are used to reconstruct a three-dimensional image.

The maximum entropy on the mean (MEM) has been rigorously used as an inverse solver. In the previous research works, the information of the neuronal sources gathered via the functional-magnetic-resonance-imaging (fMRI) technique is fed to the MEM algorithm. In Beloucha et al. [25], the authors have suggested to use diffusion-MRI (dMRI) to obtain a better visualization of the functionally connected regions of the brain.

The bidomain method, suggested in [26], iteratively converges to the optimal solution by considering the constantly active neural tissues tied with nonlinear cell model. The bidomain method relies upon the two-component model of the cardiac electrical conduction. This reaction and diffusion model considers the direct dependence and interactions between intracellular and extracellular domains [27]. Ordinary differential equations are used to define the electrical conduction model of the cell. While incorporating this approach, it is intended to obtain the brain sources that produce the scalp potential recorded by EEG calculations considering a non-static reconstruction. Unlike the conventional techniques, it does not assume any quasi-static approximation. The bidomain method has been used in the past to describe the electrical conduction model of the heart. However, recently it has been adaptively modified and used as a formidable EEG inverse solver [28–30]. It uses the lead-field matrix, in conjunction with the cell model, to define a mathematical relationship between the scalp potentials and the stimuli in the cell model. The bidomain method does not use regularization technique on the neuronal current sources like other methods; instead, it regularizes the stimuli that tend to generate these neuronal current sources.

Recently, another technique to solve the ill-posed nonlinear computed-tomography inverse problem has been proposed by synergistically combining the conventional regularization techniques with the deep-learning algorithms in Adler and Öktem [31]. The proposed scheme basically utilizes the iterative gradient-descent scheme to regularize the data. The ‘gradient’ is iteratively updated with the aid of a convolutional neural network [32]. The gradients of data discrepancy and regularizer act as the inputs to the network.

The organization of the remaining paper is as follows. Section 2 provides a comprehensive overview of the widely used conventional techniques in EEG source analysis. The explanation of the recent techniques, along with their comparative performance assessment, is presented in Section 3. Section 4 presents a detailed critical analysis of the reviewed techniques and provides recommendations for possible future enhancements. The paper is concluded in Section 5.

2. Conventional techniques

This section, the lead-field matrix and the forward problem are elaborated, followed by a brief overview of the widely used conventional EEG source localization techniques; namely, LORETA and MNE. The electrical conduction model of the brain can be described by Poisson’s equation, as given by the following equation, [33,34]. (1) $\mathrm{\nabla }.\sigma \mathrm{\nabla }\varphi =-\overline{\mathit{J}}$(1) with the boundary conditions given by the following equation. (2) $\sigma \mathrm{\nabla }\varphi .\mathit{n}=0,$(2) where $\sigma$ is denoted as conductivity, $\varphi$ is the static electrical potentials, $\overline{\mathit{J}}$ is the source distribution in conduction volume, and $\mathit{n}$ is the outward orthogonal vector. The current density’s orthogonal component is considered zero.

The pure Neumann problem, described by Equation (1), does not yield a trivial solution. However, it can be solved to deliver a unique solution by applying additional restrictions; reducing it to Laplace’s equation [35], fixing one of the scalp-electrodes as reference and setting its value zero [36,37], or employing the technique presented in [38]. As discussed earlier, the lead-field-matrix, $G$, provides a projection between the source distribution and the measured scalp potentials. The mathematical relationship for this phenomenon is given by the following equation: (3) $e=Gs+d,$(3) where $e$ is denoted as the vector of measured scalp-potentials, $s$ is denoted as the vector of active neuronal sources, and $d$ is the additive noise. The lead-field matrix is used to reconstruct a neural source image using the scalp-potentials. The expression in Equation (3) cannot be solved in straight-forward manner because the matrix $G$ is non-invertible. Hence, regularization techniques are employed to solve the NSI problem. The minimization criterion is given by the following equation: (4) $\text{min}\parallel e-Gs\parallel .$(4)

2.1. Minimum norm estimates

The MNE method is used to reconstruct the neural source image resulting from the electrical activity on the surface as well as within the brain volume. Owing to its lack of consideration of the prior restrictions, it is only capable of delivering a minimum energy solution; that is, a sub-optimal estimate of the actual source distribution in the brain [39,40]. The brain source locations thus evaluated are given by the following equation: (5) $s=(G^{T}G+\mu I{)}^{-1}{G}^{T}e,$(5) where $\mu$ is referred to as the regularization parameter and $I$ is the identity matrix. The value of $\mu$ was evaluated through the L-curve algorithm [41]. The maximum value of the parametric plot of L-curve yields the optimal regularization parameter. It is to be noted that the expression given in Equation (5) is derived from the Tinkonov regularization as the optimization of the following equation: (6) $\parallel e-Gs{\parallel }^2+\mu \parallel Is{\parallel }^2.$(6)

2.2. Low-resolution electromagnetic tomography activity

The functioning of LORETA technique is identical to MNE, except that it uses discrete Laplace operator for regularization instead of using the identity matrix [42]. As LORETA consider the connectivity given by the discrete Laplacian (${\mathrm{\Delta }}_D$), it has been rigorously used for regularization of smoothly distributed sources. The solution is given by the following equation in this case. (7) $s=(G^{T}G+\mu {{\mathrm{\Delta }}_D}^T{\mathrm{\Delta }}_D{)}^{-1}{G}^{T}e.$(7) It is to be noted that the expression given in Equation (8) is derived from the Tinkonov regularization; such that, $I={\mathrm{\Delta }}_D$, the Laplacian operator. Hence, the regularization becomes an H² semi-norm.

3. Recent advances in EEG inverse problem solution

This section briefly discusses the algorithmic advances in computing accurate solutions for NIS and inverse problem.

3.1. Differential inverse scattering methodology

A unique differential methodology that effectively utilizes the differential form of the wave as well as the first Maxwell equation (Gauss’s law for electricity) to solve the electromagnetic inverse problem in biomedical imaging has been proposed in [24]. Apart from quickly providing accurate results, the proposed differential methodology also tends to cover the limitations associated with the conventional integral based techniques. The inverse electromagnetic problem is described via three-dimensional integral equations derived from Green’s identities. To solve these three-dimensional equations, they are approximated to two-dimensional integrals in order to avoid complexities associated with Green’s function [43]. The forward and inverse integrals are given by Equations (8) and (9), respectively [24]. (8) $\text{Forward:}{E}_i\left(p\right)=E_t\left(p\right)-\frac{j{k}_b^2}{4}{\int }^L(p-q)\left(\epsilon \right(q)-{\epsilon }_b(q\left)\right)E\left(q\right)\text{d}q,$(8) (9) $\text{Inverse:}{E}_{s,j}\left(p\right)=\frac{j{k}_b^2}{4}{\int }^L(p_j-q)\left(\epsilon \right(q)-{\epsilon }_b(q\left)\right)E\left(q\right)\text{d}q,$(9) where $j$ is the location of the antenna, $p_j$ and $q$ represent the integral variable taken over the image domain, ${\epsilon }_b$ and $\epsilon$ represent the complex permittivity of the object and the background medium, $E_i$, $E_s$, and $E_t$ are the incident, scattered and total electric fields, respectively. The evaluation of Green’s function, $L$, is a complex and computationally expensive process. Hence, the following approximations are made in order to solve these equations.

• The domain is assumed to be homogeneous along one axis.

• The imaging antennas are considered as point sources.

• A background medium is used to reduce the effect of the galvanic integral term.

Due to these approximations, only a 2-D image from the dielectric properties of the domain has been reconstructed by iterative procedures utilizing either Born’s, Van Den Berg’s, Newton’s, or the gradient-based methods [44]. A three-dimensional differential form of the inverse electromagnetic problem has been proposed in Afasari et al. [24] that utilizes the differential form of the first Maxwell equation. Consequently, Green’s function is no longer needed. There is no need of assuming homogeneity along one axis. The requirement of background matching medium is also lifted off as differential form no longer relies on Green’s identities. Moreover, the differential methodology considers the accrual structural dimensions of imaging antennas instead of considering them as point sources.

First of all, an unknown area is irradiated with the aid of $E_i$ that is generated by the antennas. The field of the resulting scattered waves, $E_s$, is recorded to reconstruct the neural source images. The three-dimensional differential equations and boundary conditions governing the aforementioned wave propagation phenomenon are mathematically described by Equations (10)–(13). (10) $\text{Forward:}\left\{\begin{array}{l}\mathrm{\nabla }\times \mathrm{\nabla }\times E-k^{2}E=0,\\ BC:\left\{\begin{array}{l}n\times (E_1-E_2),\\ \epsilon E_n={\rho }_s,\end{array}\right.\end{array}\right.$(10) (11) $\text{such that},\epsilon ={\epsilon }_o{\epsilon }_r+\frac{\sigma }{j\omega },$(11) where $k$ is the wave-number, $E$ is the electric field, $BC$ are the boundary conditions, $n$ is the normal-vector on the boundary (pointing outwards), ${\rho }_s$ is the charge on the antenna surfaces, $E_n$ is the orthogonal component of the electric field on the surface, ${\epsilon }_r$ is the inhomogeneous relative permittivity, $\sigma$ is the conductivity of the medium, ${\epsilon }_o$ is the permittivity, and $\omega$ is the angular frequency. The expression presented in Equation (10) is the forward equation in the proposed differential architecture and serves to evaluate the electric field, $E$. The differential form of the first Maxwell equation, commonly referred to as Gauss’s law of electricity, is employed to derive the inverse differential equation. It is given by the following equation: (12) $\mathrm{\nabla }.\left(\epsilon E\right)=\rho ,$(12) where $\rho$ is the volume charge density. Thus, the corresponding inverse equation, and its associated boundary condition, in the proposed partial differential framework are given by the following equation: (13) $\text{Inverse}:\left\{\begin{array}{l}(\mathrm{\nabla }.E)\epsilon +E.\mathrm{\nabla }\epsilon =0,\\ BC:\epsilon =\frac{{\rho }_s}{E_n}.\end{array}\right.$(13) Together, the expressions in Equations (10) and (13) provide a complete forward-inverse set of equations to solve the NSI problem. In order to validate the proposed differential inverse equations, an array of 8 antennas (operating at 1 GHz) have been used to illuminate a real-size head model with main head tissues. The neural dielectric properties were derived using the studies provided in [45]. The head contains a sphere with radius 1 cm to emulate a bleeding with dielectric properties ${\epsilon }_r$ = 58.7 and $\sigma$ = 1.92 S/m [46]. Finally, the three-dimensional source reconstruction has been done with the aid of Finite Element Method and Born’s iteration technique [47].

3.2. MEM-based diffusion-MRI framework

A probabilistic technique to solve the NSI problem has been presented in Beloucha et al. [25]. The authors have utilized an MEM framework along with the information provided by dMRI to localize the neuronal sources. The MEM provides a good focal solution. The dMRI technique refers to the mixing of water in tissues that gives better visualization of fibre structures. It provides a better correlation of the white/gray matters with functional areas of the brain [48]. In the previous research papers using the MEM-based source localization technique, the multivariate source prelocalization is used in conjunction with the fMRI to parcellate the cortical surface [49].

The whole cortex is selected and the sources are assumed to be functionally independent and following a Gaussian probability distribution. For each subject a unique cortical surface is defined which is integrated with MEM framework. The sensors have been placed at various locations of head. Each task is repeated N times and the magnetic field (M), across the sensors, is measured. The resulting averaged signal at each sensor (m) is the arithmetic mean of M, and is given by the following equation: (14) $m=Gs+\alpha ,$(14) where $\alpha$ is the zero-mean Gaussian random variable. In MEM, the information from received signal is maximized by maximizing the mean signal entropy. The cortical surface is divided into sections by dMRI to form M different sections, {P₁, P₂, … , P_M}. Each segment has normally distributed sources and can either be inactive or active providing two possibilities for each region.

The dMRI data were collected from 11 healthy subjects using Trio-Tim system, whereas the white/gray matter was taken from Freesurfer. The M/EEG signals were recorded in magnetically shielded room. The M/EEG measurements were averaged to yield improved signal-to-noise ratio. Pre-parcellation methods were used to evaluate the connectivity profiles of respective cortical regions. Sources having identical connectivity profiles were indicated as a dMRI patch region. The number of regions inside the pre-parcels was selected such that eigenvalues of cross-correlation matrix were more than threshold value. The proposed source reconstruction technique has been compared with MNE, in the presence of additive noise. The MEM delivered superior performance as compared to MNE. The average absolute dipoles intensities (AADIs) formulation, given by Equation (15), was used to measure average dipole (distributed sources) activation over desired time window. (15) $\text{AADI}=\frac{{\sum }_{t=t_1}^{t_2}abs\left(s\right)}{\text{no. of time samples}}.$(15)

3.3. Iterative reweighted mixed-norm estimation

A source reconstruction methodology using EEG data with high spatio-temporal resolution to analyse the non-invasive electrical activity of brain has been proposed in references [22,23]. The proposed regularization technique used mixed-norm to impose a structured sparsity in space or time [50]. The spatial sparsity involves L₁ norm over blocks and Euclidian-norm-per-block to find the stationary source estimate (source at low or zero amplitude remained at the zero amplitude during whole window of interest). The non-convex optimization problems are normally solved using shrinkages [51], iterative reweighted L₁ [52], or iterative reweighted L₂ optimization [53]. The proposed irMxNE technique imposed a non-convex block-separable penalty, by combining the L_0.5 quasi-norm over blocks with the Frobenius-norm-per-block. The objective function to be minimized for the generic MxNE is given by the following equation: (16) $\hat{s}=\begin{array}{c}\text{arg min}\\ \text{s}\in {\text{R}}^{\left(SO\right)\times T}\end{array}\frac{1}{2}\parallel e-Gs{\parallel }_{\mathrm{Fro}}^2+\mu \sum _{s=1}^S\sqrt{\parallel s{\parallel }_{\mathrm{Fro}}},$(16) where $\text{Fro}$ refers to the Forbenius-norm-per-block. The source activation matrix, $\text{s}\in {\text{R}}^{\left(SO\right)\times T}$; where $S$ represents the number of source locations, $O$ represents the number of orthogonal dipoles at each source location, and $T$ is the number of sampling intervals [23]. In irMxNE, the non-convex objective function, given by Equation (17), was minimized iteratively for a sequence of distinct weighted Mx problems. The final estimate of sources is given by Equation (18). (17) ${\tilde{s}}^{\left(k\right)}=\begin{array}{c}\text{arg min}\\ \text{s}\in {\text{R}}^{\left(SO\right)\times T}\end{array}\frac{1}{2}\parallel e-G^{\left(k\right)}s{\parallel }_{\mathrm{Fro}}^2+\mu \sum _s\parallel s_l{\parallel }_{\mathrm{Fro}},$(17) (18) ${\hat{s}}^{\left(k\right)}=W^{\left(k\right)}{\tilde{s}}^{\left(k\right)},$(18) where $s_l$ is the source activation in a given block at a specific location $l$, $k$ is the iteration number, $\mu$ is the regulazation parameter, $W^{\left(k\right)}$ is diagonal matrix of the real-numbered weightages, $w^{\left(k\right)}$, that are evaluated iteratively via the following equation: (19) $w^{\left(k\right)}=\sqrt{\parallel {\hat{s}}^{(k-1)}{\parallel }_{\mathrm{Fro}}}.$(19) Then value of μ was tuneheuristically via trial-and-error to yield optimal solution.

3.4. Bidomain inverse formulation

The bidomain approach uses the knowledge of the electrical conduction that is caused by the depolarization phenomenon across the cellular membranes [54]. Although this technique was initially developed to describe the electrical conduction model of the cardiac tissue; however, recently it has been modified to comprehend the electrical conduction in the brain volume [55]. The bidomain model combines the brain’s electrical activity with the nonlinear model of the neuronal cell. It begins by solving the model of each distinct cell in the mesh, and uses the corresponding results to evaluate the electrical conduction dynamics in the brain volume. The bidomain model is defined by the expressions given in Equations (20) and (21), [56]. (20) $\mathrm{\nabla }.\left(M_{\mathrm{in}}\mathrm{\nabla }v\right)+\mathrm{\nabla }.\left(M_{\mathrm{in}}\mathrm{\nabla }{u}_{\mathrm{ex}}\right)=\gamma \left(C_t\frac{\mathrm{\partial }v}{\mathrm{\partial }t}+I_{\mathrm{ion}}(v,w)+I_{\mathrm{ap}}\right),$(20) (21) $\mathrm{\nabla }.\left(M_{\mathrm{in}}\mathrm{\nabla }v\right)=-\mathrm{\nabla }.\left(\right(M_{\mathrm{in}}+M_{\mathrm{ex}}\left)\mathrm{\nabla }{u}_{\mathrm{ex}}\right),$(21) where $v$ is the voltage across the membrane, $u_{\mathrm{ex}}$ is the potential of the extracellular region, $\gamma$ is the ratio of the membrane’s surface area to volume, $C_t$ is the capacitance across the membrane, $I_{\mathrm{ap}}$ is the applied stimulus, $I_{\mathrm{ion}}$ is the instantaneous ionic current, and $M_{\mathrm{in}}$ and $M_{\mathrm{ex}}$ are denoted as the tensors of intracellular and extracellular conductivity, respectively. The mathematical relation between the trans-membrane potential vector ($V$) and the extracellular potential vector ($U_{\mathrm{ex}}$) is obtained by solving Equation (21). The derived expression is given by the following equation: (22) $Q{U}_{\mathrm{ex}}=E{V}_{k+2},$(22) where $Q$ is a weighting matrix associated with the extracellular potential vector. The relationship between two consecutive iterations of the $V_k$ is given by the following equation: (23) $B{V}_{k+1}=A{V}_{k+2}.$(23) The extracellular potential vector is evaluated via the following equation: (24) $U_{\mathrm{ex}}=Q^{-1}E{A}^{-1}B{V}_{k+1}.$(24) The voltage distribution over the scalp was evaluated by simply taking the product of the lead-field matrix $G$ with the Dirichlet-to-Neumann operator (${A_{bb}}^{-1}$) and $U_{\mathrm{ex}}$. This product transformed the neuronal currents in the brain volume into potentials, as shown in the following equation: (25) $e=G{A_{bb}}^{-1}{U}_{\mathrm{ex}},$(25) (26) $\text{such that},e=\left(G{A_{bb}}^{-1}{Q}^{-1}E{A}^{-1}B\right)V_{k+1},$(26) (27) $\text{and rewritten as},e=P{V}_{k+1},$(27) (28) $\text{where},P=G{A_{bb}}^{-1}{Q}^{-1}E{A}^{-1}B.$(28)

The aforementioned relationship yielded the expression given by the following equation: (29) $s={A_{bb}}^{-1}{U}_{\mathrm{ex}}={A_{bb}}^{-1}{Q}^{-1}E{A}^{-1}B{V}_{k+1}.$(29)

A regularization technique was employed to obtain an optimal source location estimate. The Tikonov function, given by the following equation, was minimized to obtain an optimal estimate of the source locations. (30) $\text{mi}{\text{n}}_{\mathrm{s}}(\parallel -P\mathrm{\Delta }ts-e-P\mathrm{\Delta }t{I}_{\mathrm{ion}}{\parallel }^2+\mu \parallel I(s-s^{\mathrm{\prime }}\left){\parallel }^2\right).$(30)

The value of $\mu$, computed via the L-curve algorithm, was 0.05. The L-curve is constructed by plotting the two quantities, given in the following equation, versus each other as a curve [57]. (31) $(\parallel -P\mathrm{\Delta }ts-e-P\mathrm{\Delta }t{I}_{\mathrm{ion}}{\parallel }_2,\parallel I(s-s^{\mathrm{\prime }}\left){\parallel }_2\right).$(31) Once the electrode potentials were evaluated and the aforementioned parameters were appropriately selected, the inverse problem was solved using the following bidomain algorithm.

Algorithm: bidomain inverse formulation

Input: Normalized voltage signals recorded by sensors on the scalp, operator ${A_{bb}}^{-1}$

Output: e_n, s_n, U_e

For

Calculate the stimuli I_app with $\text{mi}{\text{n}}_{\mathrm{s}}(\parallel -P\mathrm{\Delta }ts-e-P\mathrm{\Delta }t{I}_{\mathrm{ion}}{\parallel }^2+\mu \parallel F(s-s^{\mathrm{\prime }}\left){\parallel }^2\right)$

Solve the bidomain model and calculate e_n, s_n $\mathrm{\nabla }.\left(M_{\mathrm{in}}\mathrm{\nabla }v\right)+\mathrm{\nabla }.\left(M_{\mathrm{in}}\mathrm{\nabla }{u}_{\mathrm{ex}}\right)=\gamma \left(C_t\frac{\mathrm{\partial }v}{\mathrm{\partial }t}+I_{\mathrm{ion}}(v,w)+I_{\mathrm{ap}}\right)$ $s_n={A_{bb}}^{-1}{U}_{\mathrm{ex}}$ $e_n=Gs$

Display the results

3.5. Iterative deep neural networks

As discussed earlier, this regularization technique is majorly built on the foundations of conventional gradient-descent scheme [31]. The proposed technique did not minimize the norm 2 of natural error objective function, given by the following equation: (32) $E\left(s\right)=\parallel e-Gs{\parallel }^2.$(32) Instead, it employed a variational regularization technique that not only minimized the data-discrepancy functional, given by Equation (32), but also minimized an additional regularized objective functional, given by Equation (33). (33) $\arg minE\left(s\right)=\arg min\left(D\right(s)+\mu M(s\left)\right),$(33) where $D\left(s\right)$ refers to the data-discrepancy functional and $M\left(s\right)$ is the regularization functional. The functional $D\left(s\right)$ is an affine transformation of the data log-likelihood. It represented the misfit existing in the actual and estimated data regarding the source locations. For simplicity, $D\left(s\right)$ was considered equal to $\parallel e-Gs{\parallel }^2$. The auxiliary regularization functional, $M\left(s\right)$, captured and provided a priori information regarding the estimate of potential source locations and imposed penalties on impracticable solutions [31]. A differential convex optimization technique, given by the following equation, was used to find the local minimum. (34) $s_{k+1}=s_k+\mathrm{\Delta }{s}_{k+1},$(34) (35) $\text{such that},\mathrm{\Delta }{s}_{k+1}=-\delta \mathrm{\nabla }E\left(s_k\right),$(35) where $\delta$ is the step-length between successive iterations. The gradient-descent scheme was ran iteratively until it converged to the optimal solution, yielding the minimum error functional. A convolutional neural network was used in conjunction with large sets of training data in order to learn the dynamical changes occurring in the gradient and update it at each iteration [32]. Furthermore, the gradient descent was accelerated by introducing persistent memory that allowed the algorithm to utilize the information from previous iterations of the true signals. The network had 3 layers. Each layer had 32 channels. A maximum of 10 iterations were allowed and the amount of memory allocated was 5. These hyper-parameters were heuristically selected via the trial-and-error method. Once the hyper-parameters were selected, the convolutional kernels were learned from the training data.

The qualitative comparative performance assessment of the aforementioned EEG source analysis and localization techniques is summarized in Table 1. The deductions provided in Table 1 have been carefully inferred from the numerical and graphical observations presented in the reviewed research works.

Table 1. Comparative performance assessment.

	Technique
Performance parameters	3.1	3.2	3.3	3.4	3.5
3D image reconstruction	Yes	No	No	No	No
Real-time monitoring	No	No	No	Yes	Yes
Computation	Offline	Offline	Offline	Online	Online
Priori data	No	No	No	Yes	Yes
Source regularization	Direct	Direct	Direct	Indirect stimulus	Direct
Cell-model information	None	None	None	Available	None
Nonlinearities considered	No	No	No	Yes	No
Predictive	No	No	No	Yes	Yes
Convergence speed	–	–	Fast	Slow	Fast
Approximation	Quasi-static	Quasi-static	Quasi-static	Continuum	–
Temporal resolution	Good	Good	Better	Good	Best
Spatial resolution	Better	Good	Good	Best	Good
Computation complexity	High	Low	Low	High	High

4. Discussions and future directions

This paper presents a concise review of the recent methodologies and trends in solving the NSI problem. Each technique offers a distinct set of benefits and demerits. This section presents a critical analysis of all of these techniques.

4.1. Critical analysis

The differential inverse scattering methodology improves upon the existing integral methods that could only yield two-dimensional images. The technique promises higher accuracy, three-dimensionality in image reconstruction, and low computation time. Under ideal conditions, and properly tuned parameters, the technique may very well exceed the conventional fMRI and dMRI imaging techniques. Maxwell’s equations are not exact but approximations. So for extremely strong fields and extremely short distances, there will be considerable inaccuracy. The scheme does not provide any means to compensate the effect of external noise. Its performance is prone to be degraded by the effects of exogenous additive noise sources. However, this issue can be addressed by applying the MEM algorithm to the data gathered by the proposed differential inverse scattering methodology in order to eliminate the random noises and disturbances.

The irMxNE technique, unlike its predecessors, is a fast non-convex optimization technique that solves the inverse problem by iteratively solving individual weighted MxNE problems. The irMxNE technique yields comparatively better results than the MxNE techniques. In renders improved amplitude-bias and source reconstruction. The spatio-temporal images of the sources, yielded by irMxNE, are easy to comprehend. It does not consider the number of correlated sources in its model which significantly reduces its computational complexity and convergence-time. Its iterative nature allows it to converge to the global optimum solution. However, the irMxNE is designed for offline source reconstruction and it assumes that the sources are stationary in time. Consequently, irMxNE cannot be used for real-time brain monitoring.

The bidomain formulation was also utilized as an EEG inverse solver. The reconstruction neuronal source images presented by it were temporally and spatially close to actual electrical behaviour of the sources. In contrast to other techniques, the inverse bidomain formulation incorporates the time-varying nonlinearities associated with cell model, and its effects on the conductivity in the brain. Unlike the conventional techniques that utilize the quasi-static approach, this method is capable of maintaining the continuum assumption. The bidomain formulation uses the estimated solution at a given time-instant to compute the source locations, and predicts the solution at the next time-sample. Hence, the predictive and iterative learning nature of the algorithm aids in improving the spatial and temporal resolution for source localization. Unlike irMxNE, it includes prior restrictions in the solution. The source locations estimated by the bidomain formulation relies upon several heuristically tuned parameters. This is beneficial because it increases the degree-of-freedom and the algorithm design is more flexible to attain an optimal or sub-optimal solution. All of these features enable the bidomain technique to reconstruct non-stationary sources in brain. This is extremely helpful in real-time monitoring, studying metabolic processes, and understanding the cognitive functioning of the brain. However, dependence on several tuned parameters in the bidomain may affect the solution due to the modelling errors. Optimizing all the parameters via trial-and-error is a cumbersome procedure. The bidomain formulation is a computationally complex and time-consuming process.

The iterative deep-learning based inverse solver, proposed in Adler and Öktem [31], is by far the latest and the most robust technique. The technique augments the classical regularization technique with the recent developments in deep learning. The synergistic combination of the two schemes yields the adaptability and robustness of neural networks while maintaining the optimality of the conventional regularizers. This technique also uses the prior information that is stored in the forward operator and noise model. The iterative methods in imaging are usually time-consuming, but this problem is catered in the proposed algorithm by pre-defining the maximum number of iterations. This pre-setting may reduce the flexibility of the design and restrict it from working at its full potential for a particular problem. However, these settings can be tweaked. Unlike the bidomain formulation, the technique does not rely upon an exact model of the system; instead, it iteratively learns and evolves.

4.2. Proposed future enhancements

Despite their optimality and effectiveness, there is still a lot of room for future enhancements that can be introduced in the aforementioned five techniques. The possible future directions are given as follows.

• The principle of differential inverse scattering can be extended to higher dimensions (the existing three-dimensions and time).

• The optimization rendered by the application of MEM technique on the data sets gathered by dMRI, fMRI, and differential inverse scattering can be investigated to build optimized real-time and practical microwave medical imaging techniques.

• In order to reconstruct the dynamic sources, the irMxNE technique must be retrofitted with moving-window strategies (e.g. variants of wavelet transforms), or minimization of sparsity constraints in time or frequency domain.

• The convergence rate of the bidomain inverse solver can be significantly enhanced by appropriately incorporating MxNE techniques for regularization, instead of the existing Tikhonov minimization criteria.

• The bidomain inverse solution can be utilized to classify movement intention by splitting the EEG measurements into distinct frequency bands.

• The regularization parameter (and other hyper-parameters) used in the aforementioned techniques must be optimized to their global best values. This can be done in future by investigating the meta-heuristic algorithms, evolutionary algorithms, or gradient-based optimization techniques. These techniques commonly include genetic algorithms, particle swarm optimization, bacterial foraging, artificial bee colony, gradient-descent method, etc.

• An auxiliary fuzzy inference system can be introduced alongside the convolutional neural network to enhance its capability.

• The iterative deep-learning technique should be used to reconstruct three-dimensional images using the two-dimensional training data.

• The source localization performance rendered by the integration of recurrent-wavelet neural network (instead of convolutional neural network) with classical regularizers to solve the inverse problem must be investigated.

• The scalability of using two or more of the aforementioned techniques concurrently to synergize the overall neural source image reconstruction process must be investigated.

5. Conclusion

This review paper addresses the recent advances in the field of accurately solving the NSI problem. Theoretical and mathematical backgrounds of the classical as well as the state of the art techniques have been presented. These techniques are rigorously and critically analysed. Their advantages have been properly documented. Their limitations have also been addressed along with feasible solutions to circumvent those problems. Based on the comparative analysis of these techniques, the gradient-descent method, augmented with iterative deep-learning, has proven to be quite robust and computationally efficient. There is still a lot of room for improvement in the reviewed algorithms. Therefore, a detailed list of possible future directions has also been presented to further enhance the efficiency and robustness of the recent techniques discussed in this paper.

Disclosure statement

No potential conflict of interest was reported by the authors.