A hybrid analytical-numerical algorithm for determining the neuronal current via EEG

TL;DR

This paper introduces a hybrid analytical-numerical method using radial basis functions and Tikhonov regularization to accurately reconstruct neuronal currents from EEG data, addressing the ill-posed inverse problem efficiently.

Contribution

It presents a novel hybrid algorithm combining RBF expansion and regularization for neuronal current reconstruction from EEG, improving computational efficiency and accuracy.

Findings

Reconstruction achieved RMSE of 0.1122 with synthetic data.

Method is computationally efficient and suitable for MATLAB implementation.

Effective in handling noisy EEG data with 20 dB SNR.

Abstract

In this study, the neuronal current in the brain is represented using Helmholtz decomposition. It was shown in earlier work that data obtained via electroencephalography (EEG) are affected only by the irrotational component of the current. The irrotational component is denoted by $\Psi$ and has support in the cerebrum. This inverse problem is severely ill-posed and requires that additional constraints are imposed. Here, we impose the requirement of the minimization of the $L_2$ norm of the current (energy). The function $\Psi$ is expanded in terms of inverse multiquadric radial basis functions (RBF) on a uniform Cartesian grid inside the cerebrum. The minimal energy constraint in conjunction with the RBF parametrization of $\Psi$ results in a Tikhonov regularized solution of $\Psi$. The RBF shape parameter (regularization parameter), is computed by solving a 1-D non-linear maximization problem. Reconstructions are presented using synthetic data with a signal to noise ratio (SNR) of $20$ dB. The root mean square error (RMSE) between the exact and the reconstructed $\Psi$ is RMSE=$0.1122$. The proposed reconstruction algorithm is computationally efficient and can be vectorized in MATLAB.

Tables

Table 1. Table 1: The conductivity values for the different compartments of the head model.

Domain $\Omega$	Conductivity $\sigma$ (S/m)
Cerebrum ${\Omega }_c$	0.33
CSF ${\Omega }_f$	1.0
Skull ${\Omega }_b$	0.0042
Scalp ${\Omega }_s$	0.33

Table 2. Table 2: This table shows the number of nodes and triangles in the realistic head model.

Surface meshes	Number of nodes	Number of triangles
Cerebrum: $S_c$	2562	5120
Skull $S_b$	2562	5120
Scalp $S_c$	2562	5120

Table 3. Table 3: The table shows the values of the parameters { α i , β i , 𝝉 i : 1 ≤ i ≤ 4 } conditional-set subscript 𝛼 𝑖 subscript 𝛽 𝑖 subscript 𝝉 𝑖 1 𝑖 4 \{\alpha_{i},\beta_{i},\boldsymbol{\tau}_{i}:1\leq i\leq 4\} used for generating data.

${\alpha }_i$	${\beta }_i$	${𝝉}_{x,i}$	${𝝉}_{y,i}$	${𝝉}_{z,i}$
$2$	$1e3$	$-0.0638$	$-0.0185$	$0.0546$
$-2$	$1e3$	-0.0424	0.0630	0.0324
$2$	$1e3$	0.0076	-0.0185	0.1139
$-2$	$1e3$	0.0433	-0.0366	0.0398

Equations

$$u_s(\mathbf{r})=-\frac{1}{4\pi }{\int }_{{\Omega }_c}{\nabla }_{\mathbf{\tau }}^2\Psi (\mathbf{\tau })v_s(\mathbf{r},\mathbf{\tau })dV(\mathbf{\tau }),\mathbf{r}\in S_s,$$

$${\mathbb{J}}^p(\mathbf{\tau })={\nabla }_{\mathbf{\tau }}\Psi (\mathbf{\tau })+{\nabla }_{\mathbf{\tau }}\times \mathbf{A}(\mathbf{\tau }),\mathbf{\tau }\in {\Omega }_c,$$

$$u_s(\mathbf{r})=-\frac{1}{4\pi }{\int }_{{\Omega }_c}{\nabla }_{\mathbf{\tau }}^2\Psi (\mathbf{\tau })v_s(\mathbf{r},\mathbf{\tau })dV(\mathbf{\tau }),\mathbf{r}\in S_s,$$

$$u_{s}(\boldsymbol{r})=\frac{1}{4\pi}\int_{S_{c}}\boldsymbol{n}\cdot\biggl{[}\Psi(\boldsymbol{\tau})\nabla_{\boldsymbol{\tau}}v_{s}(\boldsymbol{r},\boldsymbol{\tau})-v_{s}(\boldsymbol{r},\boldsymbol{\tau})\nabla_{\boldsymbol{\tau}}\Psi(\boldsymbol{\tau})\biggr{]}dS(\boldsymbol{\tau}),\quad\boldsymbol{r}\in S_{s}.$$

$${\nabla }_{\mathbf{\tau }}^2{v}_c(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_c,\mathbf{\tau }\in {\Omega }_c,$$

$$\displaystyle\frac{\partial}{\partial n}\biggl{[}\frac{1}{|\boldsymbol{r}-\boldsymbol{\tau}|}+v_{c}(\boldsymbol{r},\boldsymbol{\tau})\biggr{]}=\sigma_{f}\frac{\partial v_{f}(\boldsymbol{r},\boldsymbol{\tau})}{\partial n},\quad\boldsymbol{r}\in S_{c};$$

$${\nabla }_{\mathbf{\tau }}^2{v}_f(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_f,\mathbf{\tau }\in {\Omega }_c,$$

$$\displaystyle v_{f}(\boldsymbol{r},\boldsymbol{\tau})=\frac{1}{\sigma_{c}}\biggl{[}\frac{1}{|\boldsymbol{r}-\boldsymbol{\tau}|}+v_{c}(\boldsymbol{r},\boldsymbol{\tau})\biggr{]},\qquad\boldsymbol{r}\in S_{c},$$

$${\sigma }_f\frac{\partial v_f(\mathbf{r},\mathbf{\tau })}{\partial n}={\sigma }_b\frac{\partial v_b(\mathbf{r},\mathbf{\tau })}{\partial n},\mathbf{r}\in S_f,\mathbf{\tau }\in {\Omega }_c;$$

$${\nabla }_{\mathbf{\tau }}^2{v}_b(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_b,\mathbf{\tau }\in {\Omega }_c,$$

$$v_b(\mathbf{r},\mathbf{\tau })=v_f(\mathbf{r},\mathbf{\tau }),\mathbf{r}\in S_f,$$

$${\sigma }_b\frac{\partial v_b(\mathbf{r},\mathbf{\tau })}{\partial n}={\sigma }_s\frac{\partial v_s(\mathbf{r},\mathbf{\tau })}{\partial n},\mathbf{r}\in S_b,\mathbf{\tau }\in {\Omega }_c;$$

$${\nabla }_{\mathbf{\tau }}^2{v}_s(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_s,\mathbf{\tau }\in {\Omega }_c,$$

$$v_s(\mathbf{r},\mathbf{\tau })=v_b(\mathbf{r},\mathbf{\tau }),\mathbf{r}\in S_b,$$

$$\frac{\partial v_s(\mathbf{r},\mathbf{\tau })}{\partial n}=0,\mathbf{r}\in S_s.$$

$${\sigma }_c{\nabla }_{\mathbf{\tau }}^2{u}_c(\mathbf{r},\mathbf{\tau })={\nabla }_{\mathbf{\tau }}\cdot \mathbf{Q}\delta (\mathbf{r}-\mathbf{\tau }),\mathbf{r}\in {\Omega }_c,$$

$${\sigma }_c\frac{\partial u_c(\mathbf{r},\mathbf{\tau })}{\partial n}={\sigma }_f\frac{\partial u_f(\mathbf{r},\mathbf{\tau })}{\partial n},\mathbf{r}\in S_c;$$

$${\nabla }_{\mathbf{\tau }}^2{u}_f(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_f,$$

$$u_f(\mathbf{r},\mathbf{\tau })=u_c(\mathbf{r},\mathbf{\tau }),\mathbf{r}\in S_c,$$

$${\sigma }_f\frac{\partial u_f(\mathbf{r},\mathbf{\tau })}{\partial n}={\sigma }_b\frac{\partial u_b(\mathbf{r},\mathbf{\tau })}{\partial n},\mathbf{r}\in S_f;$$

$${\nabla }_{\mathbf{\tau }}^2{u}_b(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_b,$$

$$u_b(\mathbf{r},\mathbf{\tau })=u_f(\mathbf{r},\mathbf{\tau }),\mathbf{r}\in S_f,$$

$${\sigma }_b\frac{\partial u_b(\mathbf{r},\mathbf{\tau })}{\partial n}={\sigma }_s\frac{\partial u_s(\mathbf{r},\mathbf{\tau })}{\partial n},\mathbf{r}\in S_b;$$

$${\nabla }_{\mathbf{\tau }}^2{u}_s(\mathbf{r},\mathbf{\tau })=0,\mathbf{r}\in {\Omega }_s,$$

$$u_s(\mathbf{r},\mathbf{\tau })=u_b(\mathbf{r},\mathbf{\tau }),\mathbf{r}\in S_b,$$

$$\frac{\partial u_s(\mathbf{r},\mathbf{\tau })}{\partial n}=0,\mathbf{r}\in S_s.$$

$$u_j(\mathbf{r},\mathbf{\tau })=\frac{1}{4\pi }\mathbf{Q}(\mathbf{\tau })\cdot {\nabla }_{\mathbf{\tau }}{v}_j(\mathbf{r},\mathbf{\tau }),$$

$$j\in \{f,b,s\},\mathbf{r}\in {\Omega }_j,\mathbf{\tau }\in {\Omega }_c.$$

$$v_s(\mathbf{r},\mathbf{\tau })=\sum _{n=1}^{\infty }{H}_n\frac{2n+1}{4\pi }\frac{{\tau }^n}{r^{n+1}}{P}_n(\hat{\mathbf{\tau }}\cdot \hat{\mathbf{r}}),$$

$$f(\underline{\mathbf{r}}(b))-f(\underline{\mathbf{r}}(a))={\int }_C\nabla f(\underline{\mathbf{r}}(t))\cdot d\underline{\mathbf{r}}(t),$$

$$v_s(\mathbf{r},\mathbf{\tau })-v_s(\mathbf{r},0)=4\pi {\int }_0^{\tau }{u}_s(\mathbf{r},t\hat{\mathbf{\tau }},\hat{\mathbf{\tau }})dt,$$

$$E={\int }_{{\Omega }_c}|{\mathbb{J}}^p{|}^{2}dV(\mathbf{\tau }).$$

$$\displaystyle|\mathrm{\mathbb{J}}^{p}|^{2}=\biggl{(}\nabla\Psi+\nabla\times\boldsymbol{A}\biggr{)}\cdot\biggl{(}\nabla\Psi+\nabla\times\boldsymbol{A}\biggr{)}.$$

$$=|\nabla \Psi {|}^2+|\nabla \times \mathbf{A}{|}^2+2\nabla \Psi \cdot \nabla \times \mathbf{A}.$$

$$E=\int_{\Omega_{c}}\biggl{(}|\nabla\Psi|^{2}+|\nabla\times\boldsymbol{A}|^{2}+2\nabla\Psi\cdot\nabla\times\boldsymbol{A}\biggr{)}dV(\boldsymbol{\tau}).$$

$${\int }_{{\Omega }_c}\nabla \Psi (\mathbf{\tau })\cdot \nabla \times \mathbf{A}(\mathbf{\tau })dV(\mathbf{\tau })=0.$$

$$E=\int_{\Omega_{c}}\biggl{(}|\nabla\Psi|^{2}+|\nabla\times\boldsymbol{A}|^{2}\biggr{)}dV(\boldsymbol{\tau}).$$

$${\int }_{{\Omega }_c}|\nabla \Psi {|}^{2}dV(\mathbf{\tau }).$$

Full text

Keywords: Inverse Problems, electroencephalography, machine learning, surrogate modelling, artificial neural networks.

The medical significance of Electro-encephalography, EEG, is well established, see for example [1, 2, 3]. The estimation of the neuronal current from the measured electric potential (units Volts) on the surface of the head provided via the EEG data, can be formulated as a mathematical inverse problem [4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. The problem of computing the electric potential for a given head model and a given configuration of dipole sources is referred to as the forward problem [14, 15, 16, 39]. Solving the forward problem is a pre-requisite for the solution of the inverse EEG problem. If the head model is approximated as a volume conductor consisting of nested compartments with constant conductivities, then the forward problem can be formulated as a set of boundary integral equations [14, 15, 16, 17, 39, 41, 42]. OpenMEEG [16, 17] is an accurate boundary element solver (BEM) that solves these boundary integral equations. More precisely, it solves the forward problem for an arbitrary-shaped piecewise homogeneous conductor and a set of dipole sources. In this setting, the standardized low resolution brain electromagnetic tomography (sLORETA) technique [8] is known to provide accurate solutions. Other approaches which also use discrete formulations can be found in [43, 44, 45]. In contrast to the above important approaches, here the neuronal current is modelled as a continuous vectorial function. This is a more accurate representation of the underlying physics. Our study concetrates on the inverse problem of reconstructing $\Psi$ from EEG measurements. For solving this inverse problem, we use a result from [12] that expresses the electric potential on the scalp (measured via EEG) in terms of an integral over the cortical volume involving the product of the Laplacian of $\Psi$ and a certain auxiliary function denoted by $v_{s}$; the latter depends on the geometry of the various compartments and their conductivities but not on the current:

[TABLE]

where $u_{s}(\boldsymbol{r})$ denotes the electric potential (units Volts) measured on the scalp at position vector $\boldsymbol{r}$ via EEG, $\Psi(\boldsymbol{\tau})$ is the unknown scalar function to be reconstructed (irrotational component of the current), $\Omega_{c}$ denotes the volume of the cerebrum and $S_{s}$ is the surface of the scalp.

The outline of the paper is as follows: the head model considered in this study is discussed in section 1. The measurement equation, expressing the electric potential (units Volts) on the scalp in terms of $\Psi(\boldsymbol{\tau})$, and the auxiliary function $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ is reviewed in section 2. The minimum norm framework and the radial basis function (RBF) parametrization of $\Psi(\boldsymbol{\tau})$ is derived in section 3. Numerical results and reconstructions for the case of a realistic head model are presented in section 4. Finally, conclusions are discussed in section 5.

1 Head Model

The different compartments of the head model are shown in figure 1. The bounded domain $\Omega_{c}$ represent the cerebrum, which has conductivity $\sigma_{c}$. A shell $\Omega_{f}$ with conductivity $\sigma_{f}$, representing the cerebrospinal fluid, surrounds the domain $\Omega_{c}$. The cerebrospinal fluid that is surrounded by the skull is characterized by the domain $\Omega_{b}$ with conductivity $\sigma_{b}$. Finally, the skull is surrounded by the scalp that is modelled as a shell $\Omega_{s}$ with conductivity $\sigma_{s}$. Notations for the surfaces forming the boundaries of above domains are introduced in figure 1. The domain exterior to the head is denoted by $\Omega_{e}$, and it is assumed that $\Omega_{e}$ is not conductive. The permeability of all domains are equal to the permeability $\mu$ of empty space.

Table 1 presents the conductivity values of the head model as documented in [47, 48, 49].

The physiology of the head model is accurately characterized by the four layer compartments shown in figure 1, and the results derived in [12] are valid in the four layer head model. It can be observed from table 1 that $\Omega_{f}$ (CSF), has a higher conductivity than the remaining compartments, but it also has a very small thickness. The detailed analysis of [50] shows that the brain-CSF interface has a negligible effect in the analysis. For this reason, for our numerical examples we ignore $\Omega_{f}$ (CSF) and restrict our analysis to the three layer head model involving the compartments $\Omega_{c},\Omega_{b}$ and $\Omega_{s}$.

2 The Measurement Equation

Let $\mathrm{\boldsymbol{J}}^{p}(\boldsymbol{\tau})$, $\boldsymbol{\tau}\in\Omega_{c}$, denote the primary current (neuronal current) which is assumed to be supported within the cerebrum $\Omega_{c}$. Under the assumption that $\mathrm{\mathbb{J}}^{p}(\boldsymbol{\tau})$ has sufficient smoothness (continuous derivatives), we can employ the Helmholtz decomposition to express $\mathrm{\mathbb{J}}^{p}$ in terms of its irrotational and solenoidal components:

[TABLE]

where $\boldsymbol{A}(\boldsymbol{\tau})$ characterises the solenoidal component and satisfies the constraint $\nabla\cdot\boldsymbol{A}(\boldsymbol{\tau})=0$. $\Psi(\boldsymbol{\tau})$ characterises the irrotational component. Under the assumption that $\mathrm{\mathbb{J}}^{p}(\boldsymbol{\tau})$ vanishes on the cortical surface $S_{c}$, it is shown in [12] that the electric potential measured on the scalp, $S_{s}$, is given by

[TABLE]

where the auxiliary function $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ depends on geometry and conductivities but is independent of the current [12]. It was shown further in [13] that equation (3) can be expressed in the form of a surface integral :

[TABLE]

Thus, $u_{s}(\boldsymbol{r})$ is only affected by the values of $\Psi$ and of $\nabla_{\boldsymbol{\tau}}\Psi$ on the surface of $S_{c}$. It is therefore impossible to determine the $3$-dimensional nature of $\Psi$ from the knowledge of $u_{s}(\boldsymbol{r})$. Assuming that $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ can be computed, equation (3) provides a relation between $u_{s}(\boldsymbol{r})$ and $\nabla_{\boldsymbol{\tau}}^{2}\Psi(\boldsymbol{\tau})$. Similarly, assuming that both $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ and $\boldsymbol{n}\cdot\nabla_{\boldsymbol{\tau}}v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ can be computed, equation (4) provides a relation between $u_{s}(\boldsymbol{r})$, $\Psi(\boldsymbol{\tau})$ and $\nabla_{\boldsymbol{\tau}}\Psi(\boldsymbol{\tau})$.

2.1 Auxiliary Functions $v_{j}(\boldsymbol{r},\boldsymbol{\tau})$

It was shown in [12] that for a given geometry, the functions $v_{j}(\boldsymbol{r},\boldsymbol{\tau})$, $j=c,f,b,s$ are defined via the following boundary value problem:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Equations (5) - (8) are independent of the current $\mathrm{\mathbb{J}}^{p}(\boldsymbol{\tau})$ and depend only on the geometry and on the conductivities $\sigma_{c}$, $\sigma_{b}$, $\sigma_{f}$ and $\sigma_{s}$.

It is shown in [12] that the functions $v_{j}(\boldsymbol{r},\boldsymbol{\tau})$ can be related to the functions $u_{j}(\boldsymbol{r},\boldsymbol{\tau})$ with unit Volts, $\boldsymbol{r}\in\Omega_{j},~{}~{}\boldsymbol{\tau}\in\Omega_{c},~{}~{}j\in\{c,f,b,s\}$. These functions are defined in terms of a single dipole with moment $\boldsymbol{Q}(\boldsymbol{\tau})$ with unit Coulomb-meter, located at the position vector $\boldsymbol{\tau}$ via the following equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The functions $u_{j}$ and $v_{j}$ are related by the equation

[TABLE]

It was shown in [10, 11, 12] that for the the particular case of the spherical head model, $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ is given by

[TABLE]

where the coefficients $H_{n}$ depend only on the conductivities and the geometry and do not depend on the neuronal current. Equation (14) shows that the auxiliary function $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ depends on the three variables $(\tau,r,\hat{\boldsymbol{r}}\cdot\hat{\boldsymbol{\tau}})$. It is worth noting that in the particular case of the spherical head model, $v_{s}(\boldsymbol{r},\mathbb{0})=0$.

In the case of a realistic head model, OpenMEEG [16] accurately solves the boundary value problem described by equations (9)-(12) but not equations (5)-(8). However, we can employ the fundamental theorem of line integrals to approximate $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ from solutions of $u_{s}(\boldsymbol{r},\boldsymbol{\tau})=\frac{1}{4\pi}\boldsymbol{Q}\cdot\nabla_{\boldsymbol{\tau}}v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ computed by OpenMEEG.

For completeness, we recall the fundamental theorem of line integrals. Suppose that $C$ is a smooth curve given by $\underline{\mathbf{r}}(t)$, $a\leq t\leq b$ and let $f(\underline{\mathbf{r}}(t))$ be a function whose gradient $\nabla f(\underline{\mathbf{r}}(t))$ is continuous on $C$. Then,

[TABLE]

where $\underline{\mathbf{r}}(a)$ and $\underline{\mathbf{r}}(b)$ represent the initial and final points on $C$, respectively (the above theorem holds in any number of dimensions).

If we let $f=v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ and $\underline{\mathbf{r}}(t)=\hat{\boldsymbol{\tau}}t$, $0\leq t\leq\tau$ in equation (15), then we find

[TABLE]

where $v_{s}(\boldsymbol{r},\mathbb{0})$ is a constant and $4\pi u_{s}(\boldsymbol{r},\boldsymbol{\tau},\hat{\boldsymbol{\tau}})=\hat{\boldsymbol{\tau}}\cdot\nabla_{\boldsymbol{\tau}}v_{s}(\boldsymbol{r},\boldsymbol{\tau})$. It is straightforward to compute numerically the one-dimensional integral of equation (16). In the case of the spherical head model, $v_{s}(\boldsymbol{r},\mathbb{0})=0$ (see equation (14)). If the centroid of the brain mesh of a realistic head model is denoted by $\mathbb{c}$, then $v_{s}(\boldsymbol{r},\mathbb{c})$ is an additive constant and has no effect on the inversion formulas of equations (3) and (4). If we estimate $\mathbb{c}$ from the triangulated surface mesh of the brain by fitting a sphere to its nodes, then a reasonable assumption is $v_{s}(\boldsymbol{r},\mathbb{c})\approx 0$. The line integral of equation (16) has to be computed numerically.

The accuracy of reconstructing $\Psi(\boldsymbol{\tau})$ using the volume integral of equation (3) is dependent on $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$. One approach is to employ a large number of sources $\boldsymbol{Q}(\boldsymbol{\tau})$, to generate the required data set $u_{s}(\boldsymbol{r},\boldsymbol{\tau},\hat{\boldsymbol{\tau}})$ and to approximate the line integral of equation (16), using a simplified integration method. However, this approach might be computationally expensive and requires certain level of memory management. Furthermore, it is also well known that sources placed too close to the surface of cerebrum $S_{c}$, can result in large numerical errors in $u_{s}(\boldsymbol{r},\boldsymbol{\tau},\hat{\boldsymbol{\tau}})$, see [16, 17]. This in turn, will introduce errors in the upper limits of the line integral of equation (16). A computationally tractable alternative, that can mitigate these errors involves a carefully constructed regression model of $u_{s}(\boldsymbol{r},\boldsymbol{\tau},\hat{\boldsymbol{\tau}})$ using data generated via OpenMEEG. The details of this construction are discussed in section 2.2.

2.2 Regression Model of $u_{s}(\boldsymbol{r},\boldsymbol{\tau},\hat{\boldsymbol{\tau}})$ For A Realistic Head Model

Regression models constructed using data generated by PDE solvers are referred to as surrogate models [18, 19, 20, 21, 22, 23]. They are essentially machine learning models, that approximate a mapping between a set of inputs and outputs of a simulation. Suppose we are given a finite sample pairs of data, $(\mathbb{x}_{n},\mathbb{y}_{n})_{n=1}^{N}$ (the training data), where $\mathbb{x}\in\mathbb{R}^{d_{in}}$ and $\mathbb{y}\in\mathbb{R}^{d_{out}}$. These pairs of data represent the inputs and outputs of a computationally expensive simulation model $\mathbb{y}=f(\mathbb{x})$. Our goal is to construct a function $\hat{f}:\mathbb{R}^{d_{in}}\rightarrow\mathbb{R}^{d_{out}}$, which is as close as possible to the true function $f$. In this study, function $f=u_{s}(\boldsymbol{r},t\hat{\boldsymbol{\tau}},\hat{\boldsymbol{\tau}})$ on the right hand side of equation (16). A vital step in constructing an accurate surrogate model is the correct identification of the inputs to the model. In this connection we use our insight from the analysis of the spherical head model: equation (14) shows that $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ is a function of $r$, which is the radial distance from the center of the coordinate system to the sensor position of $\tau$, which is the radial distance from the center of the coordinate system to the source position vector, and of $(\hat{\boldsymbol{r}}\cdot\hat{\boldsymbol{\tau}})$, which is the cosine of the angle between the source and observation unit vectors. Thus, the necessary inputs to the surrogate model are $(r,\tau,\hat{\boldsymbol{r}}\cdot\hat{\boldsymbol{\tau}})$. An accurate learning of the function $u_{s}(\boldsymbol{r},t\hat{\boldsymbol{\tau}},\hat{\boldsymbol{\tau}})$ depends on a dense sampling of the inputs $(r,\tau,\hat{\boldsymbol{r}}\cdot\hat{\boldsymbol{\tau}})\in[r_{min},r_{max}]\times[\tau_{min},\tau_{max}]\times[-1,1]$. The position vectors of the observation points and source points from the centroid $\mathbb{c}$ of the brain are denoted by $\boldsymbol{r}$ and $\boldsymbol{\tau}$ respectively. The position vectors of the nodes of the brain with respect to $(0,0,0)$ are denoted by $\{\mathbb{x}_{j}\in\mathbb{R}^{3}:1\leq j\leq V_{b}\}$; similarly the position vectors of the nodes of the scalp are denoted by $\{\mathbb{y}_{j}\in\mathbb{R}^{3}:1\leq j\leq V_{s}\}$. Recall that OpenMEEG solves the boundary value problem described by equations (9)-(12). For any source-observation pair $(\boldsymbol{\tau},\boldsymbol{r})$, OpenMEEG employs a current dipole source with moment $\boldsymbol{Q}(\boldsymbol{\tau})$ Coulomb-meter. This solution, which has units Volts, will be denoted by $u_{s}(\boldsymbol{r},\boldsymbol{\tau},\boldsymbol{Q}(\boldsymbol{\tau}))$. The steps needed for generating the data points required for constructing the surrogate model are given below:

The above steps will result in a total of $M\times 3V_{b}$ data points. For the purpose of constructing the surrogate model it is convenient to organize the data as a matrix $\mathbb{d}\in\mathbb{R}^{3V_{b}M\times 4}$. The first three columns ($r,\tau,\hat{\boldsymbol{\tau}}\cdot\hat{\boldsymbol{r}}$) feature as inputs and the last column (electric potential) is the output.

We employ the commercial software package pSeven by datadvance [26] to construct the regression model (machine learning model). A discussion of the features of pSeven and its comparison to open source software is presented in [25].

A detailed and thorough investigation of machine-learning models (surrogate models) for the forward problem associated with EEG is outside the scope of this paper. It is work in progress and will be published elsewhere. It is noted that pSeven has a smart selection option for scanning through a set of algorithms to select the best model amongst the set. It performs for each training set a numerical optimization of the technique as well as its parameters [27, 28] by minimizing the cross-validation error, see [25]. Among the algorithms scanned by pSeven are the following: ridge regression [29], stepwise regression [30], elastic net [31], Gaussian processes [32], sparse Gaussian processes [33, 34], High Dimensional Approximation (HDA) [25, 35], and High dimensional approximation combined with Gaussian processes (HDAGP) (this technique is related to artificial neural networks and, more specifically, to the two-layer perceptron with a non-linear activation function[35]). Two desirable features of pSeven are: (i) all data manipulation is done via graphical user interface (GUI) and (ii) it can export the constructed surrogate model as a stand alone function in a number of scientific computing languages, including Matlab, C source for MEX, C source for stand alone program, C header for library, C source for library, functional mock-up interface (FMU) for Co-simulation $1.0$ and executable.

3 Minimization and a Numerical Solution of The Inverse Problem

It is clear from equations (3) and (4) that the only component of the neuronal current $\mathrm{\mathbb{J}}^{p}(\boldsymbol{\tau})$ that affects EEG data is $\Psi(\boldsymbol{\tau})$. However, as stated earlier, even this scalar function can not be computed uniquely, unless one imposes an appropriate constraint. It was shown in [5, 37] that for a spherical conductor the $L_{2}$ norm of $\mathrm{\mathbb{J}}^{p}(\boldsymbol{\tau})$ (minimal energy) yields a unique solution. Here, these results are generalized to the case of a realistic head model. Let us define the functional $E$ (energy) by

[TABLE]

Using equation (2) we find that

[TABLE]

Hence,

[TABLE]

However, it is shown in proposition $1$ and lemma $1$ of [36] that

[TABLE]

Thus,

[TABLE]

Taking into consideration that $u_{s}$ does not depend on $\boldsymbol{A}$ and that the second term of the right hand side of the above equation is positive, it follows that the minimization of (21) is equivalent to minimizing

[TABLE]

Thus, we have the following constrained minimization problem

[TABLE]

The requirement of the minimum norm of $\nabla\Psi(\boldsymbol{\tau})$ is a well documented choice of regularization, known as Tikhonov regularization [4, 29]. It forces the estimate $\Psi\in\mathbb{H}^{1}(\Omega_{c})$. The continuum form of the constrained optimization problem given by (22) does not take into account the impact of measurement noise on data $u_{s}(\boldsymbol{r})$. However, additive measurement noise is modelled in section 3.1.

3.1 Expansion of $\Psi(\boldsymbol{\tau})$ Using Radial Basis Functions

For the numerical minimization of the constrained optimization problem given in (22), radial basis functions (RBF) [46] are employed. We consider an EEG electrode cap with $M$ electrodes and discretize the domain $\Omega_{c}$ using $N$ cubic voxels. The function $\Psi(\boldsymbol{\tau})$ is expanded using inverse multiquadric RBFs [46], namely

[TABLE]

where $\{\boldsymbol{\tau}_{j}\in\Omega_{c},1\leq j\leq N\}$ are the position vectors of the centre of each voxel, and $\{\lambda_{j}:1\leq j\leq N\}$ denote the RBF coefficients. The inverse multiquadric basis function $\phi(\|\boldsymbol{\tau}-\boldsymbol{\tau}_{j}\|)$ is defined by

[TABLE]

The coefficient $\alpha$ in equation (24) is referred to as the shape parameter and needs to be estimated from appropriate data. In this study we propose a robust non-linear least squares strategy for estimating $\alpha$ from a data set $\{u(\boldsymbol{r}_{j}):1\leq j\leq M\}$ generated by a known $\Psi(\boldsymbol{\tau})$ which has support in $\Omega_{c}$. In this setting, a convenient representation of equation (23) takes the linear algebraic form

[TABLE]

where $\boldsymbol{\lambda}\in\mathbb{R}^{N}$ are the RBF coefficients. The entries of the matrix $\mathrm{A}\in\mathbb{R}^{N\times N}$ are given by

[TABLE]

If we assume that measurements are contaminated by additive noise, then the discretized form of equation (3) is given by

[TABLE]

where $\boldsymbol{u}[i]:=u(\boldsymbol{r}_{i})$, $\mathrm{V_{s}}[i,j]:=-\frac{v_{s}(\boldsymbol{r}_{i},\boldsymbol{\tau}_{j})}{4\pi N}$, $\boldsymbol{\Psi}_{L}[j]:=\nabla^{2}\Psi(\boldsymbol{\tau}_{j}),\quad j=1,...,N$ and $\mathbf{w}\in\mathbb{R}^{M}$ denotes the measurement noise vector. The only unknown in equation (27) are the discrete Laplacian values $\boldsymbol{\Psi}_{L}\in\mathbb{R}^{N}$. The action of the Laplacian operator $\mathrm{L}$ on the inverse multiquadradic basis function of (24) yields

[TABLE]

The optimization problem defined in (22) requires the computation of $\nabla\Psi(\boldsymbol{\tau})$. Partial derivative operators enable $\nabla\Psi(\boldsymbol{\tau})$ to be computed at discrete points $\{\boldsymbol{\tau}_{j}\}_{j=1}^{N}$. The operators for computing $\frac{\partial\Psi}{\partial\tau_{x}}$, $\frac{\partial\Psi}{\partial\tau_{y}}$, $\frac{\partial\Psi}{\partial\tau_{z}}$ are respectively given by

[TABLE]

[TABLE]

[TABLE]

3.2 Reconstructing $\Psi(\boldsymbol{\tau})$

The implementation of the constrained optimization problem given by (22) is well documented, see [4]. The shape parameter $\alpha$ that features in equations (26)-(31) for computing the operators $\mathrm{A}$ $\mathrm{L}$, $\mathrm{D}_{x}$, $\mathrm{Dy}$ and $\mathrm{Dz}$ plays a central role in the accuracy and stability of reconstructions. In fact, the shape parameter $\alpha$ can be interpreted as the regularization parameter [4]. Extensive numerical tests suggest that the optimal $\alpha$ is not sensitive to the choice of $\Psi(\boldsymbol{\tau})$ but it depends on the geometry of the head model, position vectors $\{\boldsymbol{\tau}_{j}\in\mathbb{R}^{3}:j=1,...,N\}$. Recall the Laplacian operator $\mathrm{L}$ given by equation (28), by letting $\boldsymbol{\Psi}_{L}=\mathrm{L}\lambda$ in equation (27), we find

[TABLE]

The model given by equation (32) is linear in $\boldsymbol{\lambda}\in\mathbb{R}^{N}$, non-linear in the shape parameter $\alpha$, and importantly it is a separable model [38]. The least squares error may be minimized with respect to $\boldsymbol{\lambda}$ and, thus it is reduced to a function of $\alpha$ only. The discretized form of the optimization problem of equation (22) is given by

[TABLE]

where $\|.\|^{2}$ is the Euclidean norm and $\mathrm{G}(\alpha)\in\mathbb{R}^{3N\times N}$ is defined as

[TABLE]

We employ the well documented weighted least squares method [38] to model the impact of the measurement noise $\mathbf{w}$ on the parameters $(\alpha,\boldsymbol{\lambda})$. The WLS method involves employing the covariance matrix of the measurement noise $\mathrm{C}_{\mathbf{w}}=\mathbb{E}[\mathbf{w}\mathbf{w}^{T}]$ in the cost function to be minimized. For convenience, we introduce the following notation :

[TABLE]

where $\mathrm{H}(\alpha)\in\mathbb{R}^{(M+3N)\times N}$, $\boldsymbol{0}\in\mathbb{R}^{N}$ and $\boldsymbol{\mathrm{b}}\in\mathbb{R}^{M+3N}$. The optimization problem of equation (33) is reduced to the following system of linear equations:

[TABLE]

The dependency of $\mathrm{H}(\alpha)$ in equation (35) on the parameter $\alpha$ is via the RBF matrices $\mathrm{L}$, $\mathrm{D}_{x}$, $\mathrm{D}_{y}$ and $\mathrm{D}_{z}$, respectively given by equations (28)-(31). The WLS error of equation (35) reads as

[TABLE]

The $\boldsymbol{\lambda}\in\mathbb{R}^{N}$ that minimizes equation (36) for a given $\alpha$ is given by

[TABLE]

The resulting WLS error is given by

[TABLE]

where $\mathbb{I}\in\mathbb{R}^{(M+3N)\times(M+3N)}$ is the identity matrix. Hence, the problem of finding the optimal value of the shape parameter $\alpha$ is reduced to the following 1-D maximization problem :

[TABLE]

The $\boldsymbol{\lambda}$ that minimize equation (36) for the optimal shape parameter $\alpha_{\textrm{opt}}$ is given by

[TABLE]

The optimal shape parameter $\alpha_{\textrm{opt}}$ is employed in equation (26) to compute the RBF matrix $\mathrm{A}$. Finally, $\boldsymbol{\hat{\Psi}}$ is given by

[TABLE]

4 Numerical Results

In the interest of reproducibility and to encourage further research in this direction, the code (surrogate model and reconstruction routines), head model (triangulated meshes) and the electrode positions are made available via Github https://github.com/parham1976. The realistic head model considered in this study, is taken from the sample data set in the OpenMEEG package. The number of nodes and triangles of the surface meshes is shown in table 2. The centroid of the brain mesh is $\mathbb{c}=[-0.0043,0.0169,0.0672]$ (this provides the center of the coordinate system for constructing the surrogate model).

4.1 Surrogate Model

In section 2, the steps for generating the data set required for constructing the surrogate model were described. Those steps will yield $3V_{b}M$ data points. In this study, we have $122$ electrodes, i.e., $M=122$ and $V_{b}=2562$ as given in table 2. So, a total of $3V_{b}M=937692$ data points are generated via OpenMEEG. The data set is split into a training and a test data set. A total of $890807$ data points were randomly selected for training and the remaining $46885$ data points are used for testing. For the problem under study and the $890807$ samples chosens as training data, the smart selection setting in pSeven package [26], chose a two layer neural network. We employ the error metrics relative mean distance measure (RDM)$\in[0,2]$, where minimum RDM$=0$, as well as the natural logarithm of the magnification factor ($\ln$(MAG)) [39, 40]. RDM is defined by

[TABLE]

where $u^{s}_{j}$, $u^{o}_{i}$ respectively, denote the electric potential (Volts) estimated by the surrogate model and the OpenMEEG package (reference). The natural logarithm of the magnification factor is given by

[TABLE]

Figure 2 depicts the comparison of the solution obtained via OpenMEEG with the surrogate model. For clarity of presentation only $1000$ points from the test data set of $46885$ have been plotted. However, the values of RDM$=0.2102$ and $\ln$[MAG]$=-0.0208$, are computed using all $46885$ samples of the test data. In Figure 2, the blue solid line with cross is the solution obtained via OpenMEEG and the red solid line with circle is the output of the surrogate model.

4.2 Reconstruction of $\Psi(\boldsymbol{\tau})$

We present numerical results using synthetic data. In order to avoid the inverse crime [4], we generate the data using a different model to that employed for reconstruction. We assume that $\Psi(\boldsymbol{\tau})$ takes the form

[TABLE]

where the values of $\{\alpha_{i},\beta_{i},\boldsymbol{\tau}_{i}:1\leq i\leq 4\}$ are given in the table 3.

The Laplacian of equation (44) is given by

[TABLE]

The steps for generating the numerical results presented in this section are given below:

The head model will have a sensor array of $M=122$ sensors. Setup a uniform Cartesian grid in the domain $\Omega_{c}$. Discretize $\Omega_{c}$ using $N=3293$ voxels. The length of the edge for each voxel is $h=0.0074$ meters. 2. 2.

Compute and store $\mathrm{V_{s}}[i,j]:=-\frac{v_{s}(\boldsymbol{r}_{i},\boldsymbol{\tau}_{j})}{4\pi N}$. The auxiliary function matrix $\mathrm{V_{s}}\in\mathbb{R}^{M\times N}$. 3. 3.

Use in equation (45) the values of $\{\alpha_{i},\beta_{i},\boldsymbol{\tau}_{i}:1\leq i\leq 4\}$ given in table 3 to compute $\{\nabla^{2}\Psi(\boldsymbol{\tau}_{j}):j=1,...,N\}$. Use in equation (27) $\mathrm{V_{s}}$ and $\boldsymbol{\Psi}_{L}[j]:=\{\nabla^{2}\Psi(\boldsymbol{\tau}_{j})\}$ to generate synthetic data, i.e., $\boldsymbol{u}=\mathrm{V_{s}}\boldsymbol{\Psi}_{L}+\mathbf{w}$. Here, we have employed additive white Gaussian noise (AWGN), with a signal to noise ratio of $20$ dB. 4. 4.

Solve the optimization problem of equation (39) to find the optimal shape parameter, $\alpha_{\textrm{opt}}$ and estimate $\boldsymbol{\hat{\Psi}}_{L}$. The optimal shape parameter is found to be $\alpha_{\textrm{opt}}=0.0169$. 5. 5.

Use $\alpha_{\textrm{opt}}$ to compute the linear operator $\mathrm{A}$ given by equation (26). The reconstruction is given by equation (41), i.e., $\boldsymbol{\hat{\Psi}}=\mathrm{A}\hat{\boldsymbol{\lambda}}$.

The matrix of auxiliary function values $\mathrm{V}_{s}\in\mathbb{R}^{122\times 3293}$. It has a rank($\mathrm{V}_{s}$)=$122$ and condition number $\kappa(\mathrm{V_{s}})=7.9804\times 10^{3}$. The matrix $\mathrm{H}(\alpha_{\textrm{opt}})\in\mathbb{R}^{10001\times 3293}$, which features in equation (33) has a rank($\mathrm{H}(\alpha_{\textrm{opt}})$)$=3293$ and a condition number $\kappa(\mathrm{H}(\alpha_{\textrm{opt}}))=8.4017\times 10^{7}$.

Figures 3(a)-3(d) depict the comparison between $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ for position vectors inside the domain $\{\boldsymbol{\tau}_{j}\in\Omega_{c}:1\leq j\leq N\}$ (internal nodes). To facilitate visualization of the function values $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ inside the domain $\Omega_{c}$, the position vectors $\boldsymbol{\tau}_{j}\in\Omega_{c}$ are divided into four non-intersecting sets. The blue solid line with cross are the function values $\boldsymbol{\Psi}$ given by equation (44) and red solid line with dots is the estimate $\boldsymbol{\hat{\Psi}}$ given by equation (41). The domain $\Omega_{c}$ has been discretized using a uniform Cartesian grid, i.e., $N=3293$ cubic voxels. Subplot (a) depicts the comparison of the function values of $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ for the set of position vectors $\{\boldsymbol{\tau}_{j}\in\Omega_{c}:\tau_{x}>0,\tau_{y}>0\}$. Subplot (b) depicts the comparison of the function values of $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ for the set of position vectors $\{\boldsymbol{\tau}_{j}\in\Omega_{c}:\tau_{x}<0,\tau_{y}>0\}$. Subplot(c) depicts the comparison of the function values of $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ for the set of position vectors $\{\boldsymbol{\tau}_{j}\in\Omega_{c}:\tau_{x}>0,\tau_{y}<0\}$. Finally, subplot (d) depicts the comparison of the function values of $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ for the set of position vectors $\{\boldsymbol{\tau}_{j}\in\Omega_{c}:\tau_{x}<0,\tau_{y}<0\}$. The root mean square error (RMSE) of $\boldsymbol{\Psi}$ and $\boldsymbol{\hat{\Psi}}$ is given by

[TABLE]

Using the analytic expression of $\Psi(\boldsymbol{\tau})$ given by equation (44) and $\boldsymbol{\hat{\Psi}}$ given by equation (41) in equation (46), we find that RMSE=$0.1122$.

Figure (4) depicts the comparison of $\boldsymbol{\Psi}$ as given by equation (44) on the surface $S_{c}$ with $\boldsymbol{\hat{\Psi}}$ given by equation (41) on the surface $S_{c}$.

The reconstructions in figure 4 are plotted using the open source fieldtrip package [51].

5 Conclusion

OpenMEEG package solves the electric potential (Volts) $u_{s}(\boldsymbol{r},\boldsymbol{\tau})=\frac{1}{4\pi}\boldsymbol{Q}(\boldsymbol{\tau})\cdot\nabla_{\boldsymbol{\tau}}v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ for a given dipole source but not $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$. To this end we propose a strategy involving 1-D line integrals to estimate $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ by constructing a surrogate model from data generated via OpenMEEG. The inputs to the surrogate model were selected based on insight provided by the analytic equation (14), i.e., $v_{s}(\boldsymbol{r},\boldsymbol{\tau})$ for the spherical head model. The surrogate model was constructed using the pSeven Datadvance surrogate modelling toolbox [26]. A total of $890807$ data points were randomly selected for training and the remaining $46885$ data points are used for testing. For this construction, the smart selection routine of Datadvance chose the high dimensional approximation [25], which is essentially a two-layer neural network. The relative mean distance measure (RDM) is $0.2102$ and the natural logarithm of the magnification factor ($\ln$(MAG)) is $-0.0208$.

The irrotational component of the neuronal current denoted by the scalar function $\Psi(\boldsymbol{\tau})$ has been parametrized using inverse multiquadric radial basis functions (RBFs) on a uniform Cartesian grid inside the cerebrum of a realistic head model. By employing synthetic data, the shape parameter denoted by $\alpha$ is estimated using a computationally efficient approach and was estimated to be $\alpha_{\textrm{opt}}=0.0169$. Furthermore, it is found that the shape parameter $\alpha$ is not sensitive to data, but depends on the configuration of the position vectors $\{\boldsymbol{\tau}_{j}\in\Omega_{c}:1\leq j\leq N\}$ of the centres of the RBF expansion. The choice of RBF and $\alpha_{\textrm{opt}}$ yields the condition number $\kappa(\mathrm{H}(\alpha_{\textrm{opt}}))=8.4017\times 10^{7}$ for the inversion matrix $\mathrm{H}(\alpha_{\textrm{opt}})$. Moreover, the regularization strategy involving minimizing the energy is in fact equivalent to Tikhonov regularization. Reconstructions are shown using synthetic data with a root mean square error (RMSE)=$0.1122$. The complete set of MATLAB files (surrogate model, inversion code and head model) as well as data set for reproducing the results in this paper are available from https://github.com/parham1976.

Acknowledgement

A.S. Fokas and P. Hashemzadeh are grateful to EPSRC and centre for mathematical imaging in healthcare (CMIH) for partial support. The grant number for this project is EP/N014588/1. C.B. Schönlieb acknowledges support from the EPSRC grant number EP/N014588/1. The authors are grateful for helpful comments, feedback and discussions with Martin Buhmann who is with the department of mathematics at the University of Giessen, as well as Jingwei Liang and Noémie Debroux who are with the Department of Applied Mathematics and Theoretical Physics, university of Cambridge.