The analytical subtraction approach for solving the forward problem in EEG

TL;DR

This paper introduces an analytical method for calculating potential integrals in EEG forward problem solutions, significantly improving accuracy while maintaining computational efficiency, especially for complex source configurations.

Contribution

The authors derive explicit formulas for potential integrals in the finite element method, reducing approximation errors and enabling efficient, accurate EEG forward modeling with complex head models.

Findings

Analytic formulas outperform numerical quadrature in accuracy.

Method is efficient for highly-eccentric sources in brain regions.

Comparable computational cost to simple numerical schemes.

Abstract

Objective: The subtraction approach is known for being a theoretically-rigorous and accurate technique for solving the forward problem in electroencephalography by means of the finite element method. One key aspect of this approach consists of computing integrals of singular kernels over the discretised domain, usually referred to as potential integrals. Several techniques have been proposed for dealing with such integrals, all of them approximating the results at the expense of reducing the accuracy of the solution. In this paper, we derive analytic formulas for the potential integrals, reducing approximation errors to a minimum. Approach: Based on volume coordinates and Gauss theorems, we obtained parametric expressions for all the element matrices needed in the formulation assuming first order basis functions defined on a tetrahedral mesh. This included solving potential integrals over triangles and tetrahedra, for which we found compact and efficient formulas. Main results: Comparison with numerical quadrature schemes allowed to test the advantages of the methodology proposed, which were found of great relevance for highly-eccentric sources, as those found in the somatosensory and visual cortices. Moreover, the availability of compact formulas allowed an efficient implementation of the technique, which resulted in similar computational cost than the simplest numerical scheme. Significance: The analytical subtraction approach is the optimal subtraction-based methodology with regard to accuracy. The computational cost is similar to that obtained with the lowest order numerical integration scheme, making it a competitive option in the field. The technique is highly relevant for improving electromagnetic source imaging results utilising individualised head models and anisotropic electric conductivity fields without imposing impractical mesh requirements.

Figures

Equations

$$a\left(u,v\right)=\int_{\Omega}\left\langle\overline{\mbox{\boldmath$\sigma$}}(\mbox{\boldmath$r$})\nabla u(\mbox{\boldmath$r$}),\nabla v(\mbox{\boldmath$r$})\right\rangle d\mbox{\boldmath$r$},$$

$$\displaystyle l(v)=-\int_{\Omega}\left\langle\overline{\mbox{\boldmath$\sigma$}}^{c}(\mbox{\boldmath$r$})\nabla u^{\infty}(\mbox{\boldmath$r$}),\nabla v(\mbox{\boldmath$r$})\right\rangle d\mbox{\boldmath$r$}$$

$$\displaystyle-\int_{\Gamma}v(\mbox{\boldmath$r$})\left\langle\overline{\mbox{\boldmath$\sigma$}}^{\infty}\nabla u^{\infty}(\mbox{\boldmath$r$}),\hat{\mbox{\boldmath$n$}}(\mbox{\boldmath$r$})\right\rangle d\mbox{\boldmath$r$}.$$

$$\mbox{\boldmath$K$}\mbox{\boldmath$u$}^{c}=\mbox{\boldmath$b$},$$

$$\displaystyle\mbox{\boldmath$K$}_{k}^{e}=\int_{\mathcal{T}_{k}}\left\langle\overline{\mbox{\boldmath$\sigma$}}(\mbox{\boldmath$r$})\nabla\mbox{\boldmath$\varphi$}(\mbox{\boldmath$r$}),\nabla\mbox{\boldmath$\varphi$}(\mbox{\boldmath$r$})\right\rangle d\mbox{\boldmath$r$},$$

$$\displaystyle\mbox{\boldmath$b$}_{k}^{v}=\frac{1}{4\pi\sigma^{\infty}}\int_{\mathcal{T}_{k}}\left\langle\overline{\mbox{\boldmath$\sigma$}}^{c}(\mbox{\boldmath$r$})\nabla f(\mbox{\boldmath$r$}),\nabla\mbox{\boldmath$\varphi$}(\mbox{\boldmath$r$})\right\rangle d\mbox{\boldmath$r$},$$

$$\displaystyle\mbox{\boldmath$b$}_{j}^{s}=\frac{1}{4\pi}\int_{T_{j}}\mbox{\boldmath$\varphi$}(\mbox{\boldmath$r$})\left\langle\nabla f(\mbox{\boldmath$r$}),\hat{\mbox{\boldmath$n$}}(\mbox{\boldmath$r$})\right\rangle d\mbox{\boldmath$r$},$$

$${\xi }_i=\frac{a_i+b_i{x}_1+c_i{x}_2+d_i{x}_3}{6V},i=1,2,3,4,$$

$$\mbox{\boldmath$\xi$}=\frac{1}{6V}\big{(}\tilde{\mbox{\boldmath$a$}}+\mbox{\boldmath$\Lambda$}^{T}\mbox{\boldmath$r$}\big{)},$$

$$\int_{\mathcal{T}}\xi_{1}^{i}\xi_{2}^{j}\xi_{3}^{k}\xi_{4}^{l}d\mbox{\boldmath$\xi$}=6V\frac{i!j!k!l!}{(i+j+k+l+3)!}.$$

$$\mbox{\boldmath$K$}^{e}=\frac{1}{6V}\int_{\mathcal{T}_{n}}\nabla_{\xi}\mbox{\boldmath$\varphi$}^{T}(\mbox{\boldmath$\xi$})\ \mbox{\boldmath$\Lambda$}^{T}\overline{\mbox{\boldmath$\sigma$}}\mbox{\boldmath$\Lambda$}\ \nabla_{\xi}\mbox{\boldmath$\varphi$}(\mbox{\boldmath$\xi$})d\mbox{\boldmath$\xi$},$$

$$\operatorname{vec}\left(\mbox{\boldmath$K$}^{e}\right)=\frac{1}{6V}\mbox{\boldmath$Q$}^{T}\left(\mbox{\boldmath$\Lambda$}\otimes\mbox{\boldmath$\Lambda$}\right)^{T}\ \operatorname{vec}(\overline{\mbox{\boldmath$\sigma$}}),$$

$$\mbox{\boldmath$Q$}=\int_{\mathcal{T}_{n}}\big{(}\nabla_{\xi}\mbox{\boldmath$\varphi$}(\mbox{\boldmath$\xi$})\otimes\nabla_{\xi}\mbox{\boldmath$\varphi$}(\mbox{\boldmath$\xi$})\big{)}d\mbox{\boldmath$\xi$},$$

$$\mbox{\boldmath$K$}^{e}=\frac{1}{36V}\mbox{\boldmath$\Lambda$}^{T}\overline{\mbox{\boldmath$\sigma$}}\mbox{\boldmath$\Lambda$}.$$

$$\mbox{\boldmath$s$}_{i}=\mbox{\boldmath$p$}_{i-1}-\mbox{\boldmath$p$}_{i+1},\ i=1,2,3,$$

$$\hat{\mbox{\boldmath$u$}}=\hat{\mbox{\boldmath$s$}}_{3},\ \ \hat{\mbox{\boldmath$v$}}=\hat{\mbox{\boldmath$w$}}\times\hat{\mbox{\boldmath$u$}},\ \ \hat{\mbox{\boldmath$w$}}=\hat{\mbox{\boldmath$s$}}_{1}\times\hat{\mbox{\boldmath$s$}}_{2},$$

$$\mbox{\boldmath$p$}_{1}=[0,0],\ \mbox{\boldmath$p$}_{2}=[s_{3},0],\ \mbox{\boldmath$p$}_{3}=[u_{3},v_{3}],$$

$$\displaystyle\mbox{\boldmath$r$}=[u,v,0]_{\mathcal{L}}=[u_{a},v_{a},0]_{\mathcal{L}_{a}},$$

$$\displaystyle\mbox{\boldmath$r$}_{0}=[u_{0},v_{0},-w_{0}]_{\mathcal{L}}=[0,0,-w_{0}]_{\mathcal{L}_{a}},$$

$$\displaystyle\mbox{\boldmath$R$}=[u-u_{0},v-v_{0},w_{0}]_{\mathcal{L}}=[u_{a},v_{a},w_{0}]_{\mathcal{L}_{a}}.$$

$$\hat{\mbox{\boldmath$m$}}_{i}=\hat{\mbox{\boldmath$s$}}_{i}\times\hat{\mbox{\boldmath$w$}},\ i=1,2,3.$$

$$t_1^0=\left[v_0(u_3-s_3)+v_3(s_3-u_0)\right]/s_1,$$

$$t_2^0=\left(u_0{v}_3-v_0{u}_3\right)/s_2,$$

$$t_3^0=v_0.$$

$$s_1^{-}=-\left[(s_3-u_0)(s_3-u_3)+v_0{v}_3\right]/s_1,$$

$$s_1^{+}=\left[(u_3-u_0)(u_3-s_3)+v_3(v_3-v_0)\right]/s_1,$$

$$s_2^{-}=-\left[s_3(s_3-u_0)+v_3(v_3-v_0)\right]/s_2,$$

$$s_2^{+}=\left(u_0{u}_3+v_0{v}_3\right)/s_2,$$

$$s_3^{-}=-u_0,s_3^{+}=s_3-u_0.$$

$$f_{2i}=\ln \left(\frac{R_i^{+}+s_i^{+}}{R_i^{-}+s_i^{-}}\right),$$

$$R_i^s=\frac{1}{{\left(R_i^0\right)}^2}\left(\frac{s_i^{+}}{R_i^{+}}-\frac{s_i^{-}}{R_i^{-}}\right),$$

$${\beta }_i={\tan }^{-1}\frac{t_i^0{s}_i^{+}}{{\left(R_i^0\right)}^2+|w_0|R_i^{+}}-{\tan }^{-1}\frac{t_i^0{s}_i^{-}}{{\left(R_i^0\right)}^2+|w_0|R_i^{-}},$$

$$\int_{S}\nabla_{S}\cdot\mbox{\boldmath$g$}\ dS=\int_{\partial S}\mbox{\boldmath$g$}\cdot\hat{\mbox{\boldmath$u$}}_{m}\ ds,\ \ \int_{S}\nabla_{S}g\ dS=\int_{\partial S}g\hat{\mbox{\boldmath$u$}}_{m}\ ds,$$

$$\int_{V}\nabla\cdot\mbox{\boldmath$g$}\ dV=\int_{\partial V}\mbox{\boldmath$g$}\cdot\hat{\mbox{\boldmath$u$}}_{n}\ dS,\ \ \int_{V}\nabla g\ dS=\int_{\partial V}g\hat{\mbox{\boldmath$u$}}_{n}\ dS.$$

$$\nabla_{S}\cdot\left(h(R)\frac{\mbox{\boldmath$t$}}{t^{2}}\right)=\frac{\mbox{\boldmath$t$}}{t^{2}}\cdot\nabla_{S}h(R).$$

$$\displaystyle\frac{1}{R^{3}}=-\nabla_{S}\cdot\left(\frac{1}{R}\frac{\mbox{\boldmath$t$}}{t^{2}}\right),$$

$$\displaystyle\frac{3}{R^{5}}=-\nabla_{S}\cdot\left(\frac{1}{R^{3}}\frac{\mbox{\boldmath$t$}}{t^{2}}\right),$$

Full text

Approach: Based on volume coordinates and Gauss theorems, we obtained parametric expressions for all the element matrices needed in the formulation assuming first order basis functions defined on a tetrahedral mesh. This included solving potential integrals over triangles and tetrahedra, for which we found compact and efficient formulas.

Main results: Comparison with numerical quadrature schemes allowed to test the advantages of the methodology proposed, which were found of great relevance for highly-eccentric sources, as those found in the somatosensory and visual cortices. Moreover, the availability of compact formulas allowed an efficient implementation of the technique, which resulted in similar computational cost than the simplest numerical scheme.

Significance: The analytical subtraction approach is the optimal subtraction-based methodology with regard to accuracy. The computational cost is similar to that obtained with the lowest order numerical integration scheme, making it a competitive option in the field. The technique is highly relevant for improving electromagnetic source imaging results utilising individualised head models and anisotropic electric conductivity fields without imposing impractical mesh requirements.

1 Introduction

The forward problem in electroencephalography (EEG-FP) is a cornerstone of quantitative brain activity characterisation based on electromagnetic recordings. It consists in computing the electric potential function in the head generated by a set of known current generators representing time-varying charge distributions in the cortex at the mesoscopic level. Mathematically, the EEG-FP is described by a Poisson-like equation (subject to a Neumann boundary condition) depending on a particular combination of source and conductor model parameters [1]. It is widely accepted in the literature to consider dipoles as the generators of EEG signals [2, 3, 4]. Under such assumption, analytical solutions to the EEG-FP can be obtained for multilayered spherical head models with piecewise anisotropic electrical conductivity [5]. However, it is acknowledged that the adoption of individualised head representations confer noticeable advantages over spherical models in the analysis EEG (either invasive [6] or non-invasive [7, 8, 9, 10]), leading to the use of numerical techniques.

Several numerical methodologies were presented for solving the EEG-FP, from which the finite element method (FEM) stands out of the rest. The reason is based on its flexibility for incorporating arbitrary geometries and heterogeneous and anisotropic electrical conductivity fields, which were shown to impact the current flow through biological tissues [11, 6, 12]. One key aspect in the solution of the EEG-FP by means of the FEM (as well as with any other numerical method) is related to the singularity introduced by the dipolar source model. Within few existing options for handling such singularities, the subtraction approach is the only that allows the use of standard FE basis functions for simulating dipoles in arbitrary locations while guaranteeing the existence and uniqueness of the solution [13, 1, 14]. Subtraction-based techniques were shown to provide highly accurate solutions to the EEG-FP, improving those obtained by the partial integration method, and comparable to using the Saint Venant’s principle [15, 16, 17]. Solving the EEG-FP by means of the subtraction FEM requires the solution of surface and volume integrals with singular integrands, usually referred to in electromagnetism as potential integrals [18]. So far, researchers have managed to work with approximations of these integrals obtained either projecting the singular function [14] or its gradient [13] on the FE space, or utilising numerical quadrature schemes [19, 1]. However, such approximations are known for introducing unwanted errors into the solution, requiring high resolution meshes and/or prohibitive computational resources to achieve accurate results.

In this paper, we derive analytical expressions for all the integrals needed in the subtraction FE formulation of the EEG-FP utilising linear basis functions defined on a tetrahedral mesh. This is done by exploiting the volume coordinate system and both surface and volume Gauss theorems to reduce integrals to 1D, where they become analytically solvable. Simple and compact expressions are obtained for all the integrals involved in the computation of the stiffness matrix and source vectors, avoiding any error other than that given by the domain discretisation itself. This makes the proposed method, coined analytical subtraction (AS) approach, the optimal with regard to accuracy. We illustrate the advantages of the AS technique over the use of numerical quadrature formulas at the element level, and demonstrate how these differences impact in the overall solution of the EEG-FP. Additionally, we show that the expressions allow to reduce the computational requirements compared to the use of a numerical quadrature scheme, making it a competitive option.

The rest of the paper is organised as follows: in Section 2 we review the subtraction formulation of the EEG-FP, and present the integrals needed in the FE discretisation. The solution of these integrals require different tools, and are therefore treated separately. Analytical expressions for the element stiffness matrix are derived in Section 3, whereas the corresponding to the source vector are obtained in Section 4. Local and global experiments are described in Section 5, and the corresponding results presented in Section 6. Finally, we discuss the relevance of the methodology in Section 7, and summarise the conclusions in Section 8.

2 EEG forward problem

2.1 Differential formulation

The EEG-FP consists of finding the electric potential function u(\mbox{\boldmathr}) due to a current source with density s(\mbox{\boldmathr}) defined over the domain $\Omega$ (i.e. the head), with boundary $\Gamma$. Let \overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr}) be the rank-2 conductivity tensor field within $\Omega$, and \hat{\mbox{\boldmathn}}(\mbox{\boldmathr}) the unitary vector normal to $\Gamma$. Under generally accepted assumptions (as the quasistatic and the point electrode model approximations), the EEG-FP reduces to find u(\mbox{\boldmathr}) satisfying \nabla\cdot\big{(}\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr})\nabla u(\mbox{\boldmathr})\big{)}=-s(\mbox{\boldmathr}) (\mbox{\boldmathr}\in\Omega) together with the boundary condition \left\langle\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr})\nabla u(\mbox{\boldmathr}),\hat{\mbox{\boldmathn}}(\mbox{\boldmathr})\right\rangle=0 (\mbox{\boldmathr}\in\Gamma) [5, 1]. In the case of assuming a dipolar source located in \mbox{\boldmathr}_{0} with dipolar moment $q$, s(\mbox{\boldmathr})=-\left\langle\mbox{\boldmathq},\nabla\delta(\mbox{\boldmathr}-\mbox{\boldmathr}_{0})\right\rangle.

The subtraction method consists of avoiding the singularity in s(\mbox{\boldmathr}) by separating $\Omega$ into two subsets, one surrounding the source, namely $\Omega^{\infty}$, with homogeneous conductivity \overline{\mbox{\boldmath\sigma}}^{\infty}=\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr}_{0}), and the other being $\Omega^{c}=\Omega\backslash\Omega^{\infty}$ (the complement of $\Omega^{\infty}$ in $\Omega$) with electrical conductivity \overline{\mbox{\boldmath\sigma}}^{c}(\mbox{\boldmathr})=\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr})-\overline{\mbox{\boldmath\sigma}}^{\infty}. This allows to express the electric potential as the sum of two terms, u(\mbox{\boldmathr})=u^{c}(\mbox{\boldmathr})+u^{\infty}(\mbox{\boldmathr}), where u^{\infty}(\mbox{\boldmathr}) is the singularity potential generated by a source in an unbounded homogeneous conductor with conductivity \overline{\mbox{\boldmath\sigma}}^{\infty} (for which analytical expressions exist), and u^{c}(\mbox{\boldmathr}) is the correction potential satisfying \nabla\cdot\big{(}\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr})\nabla u^{c}(\mbox{\boldmathr})\big{)}=-\nabla\cdot\big{(}\overline{\mbox{\boldmath\sigma}}^{c}(\mbox{\boldmathr})\nabla u^{\infty}(\mbox{\boldmathr})\big{)} (\mbox{\boldmathr}\in\Omega) and subject to \left\langle\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr})\nabla u^{c}(\mbox{\boldmathr}),\hat{\mbox{\boldmathn}}(\mbox{\boldmathr})\right\rangle=-\left\langle\overline{\mbox{\boldmath\sigma}}(\mbox{\boldmathr})\nabla u^{\infty}(\mbox{\boldmathr}),\hat{\mbox{\boldmathn}}(\mbox{\boldmathr})\right\rangle (\mbox{\boldmathr}\in\Gamma) [14, 19, 1]. The problem then results in finding the correction potential by means of a numerical method, after which the singularity potential is added.

2.2 FE discretisation

The FE formulation of the EEG-FP relies on the variational form of the subtraction version. This can be obtained by multiplying the corresponding differential equation by a test function $v$ belonging to a suitable space $H$, and then integrating over $\Omega$. After utilising the divergence theorem and the boundary condition [14, 19], the variational formulation results in finding u^{c}(\mbox{\boldmathr})\in H such that, for all v(\mbox{\boldmathr})\in H, satisfies $a\left(u^{c},v\right)=l(v)$, where $a:H\times H\rightarrow\mathbb{R}$ is the bilinear form given by

[TABLE]

and $l:H\rightarrow\mathbb{R}$ is the linear form defined as

[TABLE]

To proceed with the FEM, a discretisation $\mathcal{T}$ of $\Omega$ is needed. This discretisation is composed by a set of nodes \mbox{\boldmathp}_{i} ($i=1,\dots,N$) and elements $\mathcal{T}_{j}$ ($j=1,\dots,N_{e}$) defined upon these nodes. This tessellation allows to construct a discretised FE space $V_{N}\subset H$ where to compute the numerical solution. More explicitly, we choose V_{N}=span\{\varphi_{i}(\mbox{\boldmathr}):i=1,\dots,N\}, with \varphi_{i}(\mbox{\boldmathr}) being piecewise functions such that \varphi_{i}(\mbox{\boldmathp}_{j})=\delta_{ij} [14]. Then, we look for \widetilde{u}^{c}(\mbox{\boldmathr})\in V_{N} (an approximation of u^{c}(\mbox{\boldmathr})\in H) satisfying, for all v(\mbox{\boldmathr})\in V_{N}, $a(\widetilde{u}^{c},v)=l(v)$. This leads to solve the system of linear equations

[TABLE]

where \mbox{\boldmathK}\in\mathbb{R}^{N\times N} is the stiffness matrix defined by K_{ij}=a(\varphi_{i}(\mbox{\boldmathr}),\varphi_{j}(\mbox{\boldmathr})), \mbox{\boldmathb}\in\mathbb{R}^{N} is the source vector with elements b_{i}=l(\varphi_{i}(\mbox{\boldmathr})), and \mbox{\boldmathu}^{c}\in\mathbb{R}^{N} is the vector containing the numerical approximation of the correction potential on the mesh nodes (i.e. u^{c}(\mbox{\boldmathr})\approx\widetilde{u}^{c}(\mbox{\boldmathr})=\sum_{i=1}^{N}\varphi_{i}(\mbox{\boldmathr})u^{c}_{i}). The most common scenario available in the literature (and the one adopted in this work) is to utilise linear basis functions defined on a tetrahedral discretisation. In this case, the arrays $K$ and $b$ needed in (3) are obtained by assembling the corresponding element matrices, i.e. those found by integrating over individual tetrahedrons $\mathcal{T}_{j}$ ($j=1,\dots,N_{e}$) or surface triangles $T_{k}$ ($k=1,\dots,N_{s}$). For simplicity, we split the source vector into two terms representing the surface and volume integrals, i.e. \mbox{\boldmathb}=-(\mbox{\boldmathb}_{s}+\mbox{\boldmathb}_{v}), allowing their individual analysis.

Without loss of generality, we assume that $\Omega^{\infty}$ is an isotropic medium, i.e. \overline{\mbox{\boldmath\sigma}}^{\infty}=\sigma^{\infty}\mbox{\boldmathI}_{3}, with \mbox{\boldmathI}_{n} being the $n\times n$ identity matrix. Moreover, we consider dipolar source models to represent current generators, as is standard in the literature [3, 2]. In this case, u^{\infty}(\mbox{\boldmathr})=(4\pi\sigma^{\infty})^{-1}f(\mbox{\boldmathr}), with f(\mbox{\boldmathr})=\mbox{\boldmathq}\cdot\mbox{\boldmathR}R^{-3}, \mbox{\boldmathR}=\mbox{\boldmathr}-\mbox{\boldmathr}_{0}, and R=|\mbox{\boldmathR}|. Under these assumptions, the element matrices are

[TABLE]

where $k=1,\dots,N_{e}$, $j=1,\dots,N_{s}$, and \mbox{\boldmath\varphi}(\mbox{\boldmathr}) is the linear interpolation vector function with elements \varphi_{i}(\mbox{\boldmathr}), $i=1,\dots,d_{i}$ [$d_{i}=4$ for (4) and (5); $d_{i}=3$ for (6); see (19)]. In the following sections, we derive analytical expressions for (4)–(6) based on the application of linear transforms and Gauss theorems.

3 Stiffness matrix

The computation of the element stiffness matrix is easily performed by employing the volume coordinate system [20, 21]. Let $V$ be the volume of an arbitrary tetrahedron $\mathcal{T}$ with nodes \mbox{\boldmathp}_{i} ($i=1,\dots,4$). For any point \mbox{\boldmathr}\in\mathcal{T}, we define the volume coordinates $\xi_{i}=V_{i}/V$ ($i=1,\dots,4$), where $V_{i}$ is the volume of the tetrahedron with nodes $r$ and \mbox{\boldmathp}_{j}, with $j\neq i$ ($j=1,\dots,4$). These coordinates are given by

[TABLE]

where \mbox{\boldmathr}=[x_{1},x_{2},x_{3}]^{T} and $a_{i}$, $b_{i}$, $c_{i}$, and $d_{i}$ are real coefficients depending on \mbox{\boldmathp}_{i} [20]. In matrix form,

[TABLE]

where \tilde{\mbox{\boldmatha}}, \tilde{\mbox{\boldmathb}}, \tilde{\mbox{\boldmathc}}, and \tilde{\mbox{\boldmathd}} are the $4\times 1$ vectors with elements $a_{i}$, $b_{i}$, $c_{i}$, and $d_{i}$, respectively ($i=1,\dots,4$), \mbox{\boldmath\Lambda}=[\tilde{\mbox{\boldmathb}},\tilde{\mbox{\boldmathc}},\tilde{\mbox{\boldmathd}}]^{T}, and \mbox{\boldmath\xi}=[\xi_{1},\xi_{2},\xi_{3},\xi_{4}]^{T}. It is easy to note that (7) transforms $\mathcal{T}$ into a normalised tetrahedron $\mathcal{T}_{n}$ in the $\xi$ coordinate system [20].

The introduction of the volume coordinate system has two main purposes: first, it allows a simple description of polynomial basis functions of any order, reducing to \mbox{\boldmath\varphi}(\mbox{\boldmath\xi})=\mbox{\boldmath\xi} for first order FEM [21]; and second, it simplifies the solution of integrals of the form

[TABLE]

After transforming (4) to the volume coordinate system, and noting that the gradient of any scalar function h(\mbox{\boldmathr}) can be expressed in the volumetric coordinate system as \nabla h(\mbox{\boldmathr})=(6V)^{-1}\mbox{\boldmath\Lambda}\nabla_{\xi}h(\mbox{\boldmath\xi}) [21], we obtain

[TABLE]

where \nabla_{\xi}\mbox{\boldmath\varphi}(\mbox{\boldmath\xi})=\big{[}\nabla_{\xi}\varphi_{1}(\mbox{\boldmath\xi}),\dots,\nabla_{\xi}\varphi_{N}(\mbox{\boldmath\xi})\big{]} and $\nabla_{\xi}$ is the gradient operator in the volume coordinate space. To further simplify (9), we apply the $vec$ operator to \mbox{\boldmathK}^{e}, i.e. the linear operator that transforms the element matrix in a column vector obtained by stacking the columns of \mbox{\boldmathK}^{e} on top of one another. Using the identity \operatorname{vec}\left(\mbox{\boldmathA}\mbox{\boldmathB}\mbox{\boldmathC}\right)=\big{(}\mbox{\boldmathC}^{T}\otimes\mbox{\boldmathA}\big{)}\operatorname{vec}(\mbox{\boldmathB}) [22], we get

[TABLE]

where

[TABLE]

is a matrix with constant coefficients, and $\otimes$ is the Kronecker product. Once the basis functions are chosen, $Q$ is easily computed using (8). In the particular case of using first order basis functions, \nabla_{\xi}\mbox{\boldmath\varphi}(\mbox{\boldmath\xi})=\mbox{\boldmathI}_{4}, and consequently \mbox{\boldmathQ}=\mbox{\boldmathI}_{16}/6. Therefore, the element matrix considering first order basis functions turns out to be

[TABLE]

4 Source vector

The non-linearities introduced by f(\mbox{\boldmathr}) impose new difficulties in the computation of the element source vectors. To overcome them, we employ Gauss theorems defined on local reference frames to reduce their dimensionality to 1D. In this Section, we first review and extend such framework [18, 23], for then utilising it in the computation of (5) and (6).

4.1 Local coordinate system

Let $\mathcal{G}\equiv[x,y,z]$ be the global Cartesian frame in which our model is inscribed (see Fig. 1). We consider a triangle $T$ formed by the nodes \mbox{\boldmathp}_{1}, \mbox{\boldmathp}_{2}, and \mbox{\boldmathp}_{3}, such that $\partial T_{i}$ represents its $ith$ side, i.e. opposite to \mbox{\boldmathp}_{i} ($i=1,2,3$). For this triangle, we define the vectors

[TABLE]

with \mbox{\boldmathp}_{4}\equiv\mbox{\boldmathp}_{1} and \mbox{\boldmathp}_{-1}\equiv\mbox{\boldmathp}_{3}. These vectors allow to define a local frame $\mathcal{L}\equiv[u,v,w]$, with origin in \mbox{\boldmathp}_{1}, such that $T$ belongs to the plane $w=0$. The unitary vectors for the local frame are chosen as

[TABLE]

where \hat{\mbox{\boldmaths}}_{i}=\mbox{\boldmaths}_{i}/s_{i}, with s_{i}=|\mbox{\boldmaths}_{i}|, $i=1,2,3$. In this new frame, any point \mbox{\boldmathr}\in T will have coordinates $[u,v,0]_{\mathcal{L}}\equiv[u,v]$. In particular,

[TABLE]

with u_{3}=-\mbox{\boldmaths}_{2}\cdot\hat{\mbox{\boldmathu}} and v_{3}=-\mbox{\boldmaths}_{2}\cdot\hat{\mbox{\boldmathv}}.

We are interested in solving integrals involving the vector $R$ and/or its magnitude $R$. In the local frame, the source position can be written as \mbox{\boldmathr}_{0}=\mbox{\boldmath\rho}-w_{0}\hat{\mbox{\boldmathw}}, where \mbox{\boldmath\rho}=[u_{0},v_{0},0]_{\mathcal{L}} is the projection of \mbox{\boldmathr}_{0} onto the plane $w=0$. To perform such integrals, it is convenient to translate $\mathcal{L}$ so that the origin is $\rho$. This new frame, called $\mathcal{L}_{a}\equiv[u_{a},v_{a},w]$, has the same unitary vectors as $\mathcal{L}$. Noting that $u=u_{a}+u_{0}$ and $v=v_{a}+v_{0}$, it can be easily seen that

[TABLE]

The use of $\mathcal{L}_{a}$ facilitates the parameterisation of any point in $\partial T$, which will be required for calculating the line integrals after the application of Gauss theorems. Let \mbox{\boldmatht}_{i}^{0} be the shortest vector going from $\rho$ to $\partial T_{i}$. These vectors can be written as \mbox{\boldmatht}_{i}^{0}=t_{i}^{0}\hat{\mbox{\boldmathm}}_{i}, with

[TABLE]

The coefficients $t_{i}^{0}$ can be obtained from the intersection between the lines (\mbox{\boldmath\rho},\hat{\mbox{\boldmathm}}_{i}) and (\mbox{\boldmathp}_{i-1},\hat{\mbox{\boldmaths}}_{i}), resulting in [23]

[TABLE]

Then, every point \mbox{\boldmaths}\in\partial T_{i} can be parameterised as \mbox{\boldmaths}(s)=\mbox{\boldmatht}_{i}^{0}+s\hat{\mbox{\boldmaths}}_{i}, where $s\in[s_{i}^{-},s_{i}^{+}]$ satisfies \mbox{\boldmaths}(s_{i}^{-})=\mbox{\boldmathp}_{i-1} and \mbox{\boldmaths}(s_{i}^{+})=\mbox{\boldmathp}_{i+1}. After some algebra, we find [23]

[TABLE]

In the following, $R_{i}^{0}=\sqrt{(t_{i}^{0})^{2}+w_{0}^{2}}$ is the distance between \mbox{\boldmathr}_{0} and $\partial T_{i}$, and $R_{i}^{\pm}=\sqrt{(s_{i}^{\pm})^{2}+(t_{i}^{0})^{2}+w_{0}^{2}}$ is the distance between \mbox{\boldmathr}_{0} and \mbox{\boldmaths}(s_{i}^{\pm}), $i=1,2,3$. Also, to express results concisely, we define the following functions (some of them introduced in [18, 23])

[TABLE]

with the arctangent function evaluated on its principal branch.

4.2 Gauss theorems

We are interested in the application of Gauss theorems for reducing the dimensionality of 2D and 3D integrals to 1D, where analytical expressions can be easily found. Let $S$ be a planar surface, and \hat{\mbox{\boldmathu}}_{m} the unitary vector normal to its boundary $\partial S$ and belonging to the same plane (pointing outwards). Then, for any continuously differentiable functions $g:S\rightarrow\mathbb{R}$ and \mbox{\boldmathg}:S\rightarrow\mathbb{R}^{3}, the surface Gauss theorems state that [24]

[TABLE]

where \nabla_{S}=\hat{\mbox{\boldmathu}}\partial/\partial u+\hat{\mbox{\boldmathv}}\partial/\partial v is the surface vector differential operator. Similarly, let $V$ be a volume and \hat{\mbox{\boldmathu}}_{n} the unitary vector normal to $\partial V$, the boundary of $V$ (pointing outwards). Then, for any functions $g:V\rightarrow\mathbb{R}$ and \mbox{\boldmathg}:V\rightarrow\mathbb{R}^{3}, the volume Gauss theorems are given by [24]

[TABLE]

In the present study, since we are assuming tetrahedral meshes, the application of (13) to volume integrals will result in transformed integrals over their triangular faces. Therefore, special emphasis should be put in expressing the integrands in a proper differential form, so that (12) can be utilised. To this end, the following proposition (shown in Appendix 1) will be of practical use.

Proposition 1

Let $h:S\rightarrow\mathbb{R}$ be a function defined for any vector $t$ belonging to the planar surface $S$. Then, the following relation holds

[TABLE]

Proposition 1 allows to express integrands in the form required to apply Gauss theorems on a planar surface. In particular, if we choose $h(R)=R^{-1}$ and $h(R)=R^{-3}$, we get

[TABLE]

respectively. In the special case of choosing $t$ as \hat{\mbox{\boldmathu}} or \hat{\mbox{\boldmathv}}, these expressions turn out to be

[TABLE]

where the braces indicate that the relation is valid for either the upper or lower elements in them. Finally, utilising simple differentiation rules, we obtain

[TABLE]

4.3 Surface integral

Now, we take advantage of the results from the previous subsections to compute the surface element vector \mbox{\boldmathb}^{s}. In this case, the linear interpolation vector is \mbox{\boldmath\varphi}(\mbox{\boldmathr})=[\varphi_{1}(\mbox{\boldmathr}),\varphi_{2}(\mbox{\boldmathr}),\varphi_{3}(\mbox{\boldmathr})]^{T}, with \varphi_{i}(\mbox{\boldmathp}_{j})=\delta_{ij}, for every $i,j=1,2,3$. In the $\mathcal{L}$ frame, $\varphi$ takes the form [25]

[TABLE]

which clearly complies the interpolation conditions (11). Replacing (19) in (6) leads to

[TABLE]

where

[TABLE]

Then, we need to find $I_{0}$, $I_{u}^{a}$, and $I_{v}^{a}$ for any arbitrary dipolar source. To do so, the following proposition (shown in Appendix 2) is of great help.

Proposition 2

For any triangle $T$ and dipolar source located in \mbox{\boldmathr}_{0} with moment $q$, the following relation holds

[TABLE]

Proposition 2 is used to rewrite (21) – (22) as

[TABLE]

which allows to exploit the fact that, in $\mathcal{L}_{a}$, \hat{\mbox{\boldmathn}}=\hat{\mbox{\boldmathw}}, resulting in \hat{\mbox{\boldmathn}}\cdot\mbox{\boldmathR}=w_{0}.

To solve (23) in the $\mathcal{L}_{a}$ frame, we first note that, for any function $g$ defined in $\mathcal{L}_{a}$, \nabla g=\nabla_{S}g+\partial g/\partial w_{0}\hat{\mbox{\boldmathw}}. Applying this property to the integrand of (23), we obtain

[TABLE]

Eq. (25) allows us to split $I_{0}$ in three terms, which are treated separately. The first term can be easily solved by applying the second identity from (12), resulting in

[TABLE]

The computation of the second and third terms is done using (14) and (15), after which the first relation in (12) is applied. However, it should be noted that the utilisation of (14) and (15) introduce a singularity in $\rho$ (since $t=0$ for $\rho$ in $\mathcal{L}_{a}$), which prevents the application of Gauss theorems straightforwardly. This issue can be solved by splitting the domain in two parts, one being $T_{\epsilon}$, a circle with radius $\epsilon$ centred in $\rho$, and the other its complement in $T$, i.e. $T-T_{\epsilon}$. Then, the corresponding integrals are computed separately, after which the limit for $\epsilon\rightarrow 0$ is performed. This procedure was followed by Graglia [23], who obtained the following expression for the second term

[TABLE]

where \mbox{\boldmath\beta}=\beta_{1}+\beta_{2}+\beta_{3}. In Appendix 3 we show that the integral corresponding to third term is given by

[TABLE]

Finally, the proper combination of (26) – (28) yields

[TABLE]

The integrals $I^{a}_{u}$ and $I^{a}_{v}$ are obtained following a similar procedure. In Appendix 4 we show that

[TABLE]

where $R_{i}^{d}=(R_{i}^{-})^{-1}-(R_{i}^{+})^{-1}$ and $\operatorname{sgn}$ is the sign function.

4.4 Volume integral

The computation of \mbox{\boldmathb}^{v} turns out to be much simpler than the corresponding to the surface element vector. Utilising (13) and the fact that \nabla\mbox{\boldmath\varphi}=(6V)^{-1}\mbox{\boldmath\Lambda} for first order basis functions (see Section 3), (5) reduces to

[TABLE]

where $T^{i}$ is the $ith$ triangular face of the tetrahedron under consideration. Noting that, for a dipolar source, f(\mbox{\boldmathr})=-\mbox{\boldmathq}\cdot\nabla R^{-1}, we get

[TABLE]

where the last equality follows from [23, eq. (34)] (and noting that, in this paper, \mbox{\boldmathw}_{0} is opposite to its analogous in [23]) and the superscripts refer to the triangular face under consideration.

5 Experiments

5.1 Local errors

We first studied the performance of the AS method for computing the element vectors \mbox{\boldmathb}^{s} and \mbox{\boldmathb}^{v}. This was done by calculating the analytical solution for a given element and source of electrical activity, and comparing the results with those obtained with the Gauss-Jacobi quadrature scheme, known as the full subtraction (FS) approach [19, 1]. More specifically, we considered a source in \mbox{\boldmathr}_{0}=[d,0,0] with moment \mbox{\boldmathq}=[10,0,0]nAm and calculated the element vectors for equilateral simplices located in the origin with side length $a$. In the case of \mbox{\boldmathb}^{s}, we used a triangle with vertices \mbox{\boldmathp}_{1}=[0,-a/2,-a\sqrt{3}/6], \mbox{\boldmathp}_{2}=[0,a/2,-a\sqrt{3}/6], and \mbox{\boldmathp}_{3}=[0,0,a\sqrt{3}/3]; in the case of \mbox{\boldmathb}^{v}, we used a tetrahedron with nodes \mbox{\boldmathp}_{1}=[0,0,a\sqrt{3}/3], \mbox{\boldmathp}_{2}=[0,a/2,-a\sqrt{3}/6], \mbox{\boldmathp}_{3}=[0,-a/2,-a\sqrt{3}/6], and \mbox{\boldmathp}_{4}=[-a\sqrt{3}/3,0,0]. Then, we computed the relative error (RE) between element vectors, defined as RE_{e}=\|\mbox{\boldmathb}_{n}-\mbox{\boldmathb}_{a}\|/\|\mbox{\boldmathb}_{a}\|, with \mbox{\boldmathb}_{n} and \mbox{\boldmathb}_{a} being the numerical and analytical solutions, respectively. We repeated such experiments for different combinations of $d$ and $a$.

5.2 Spherical head model

In a second experiment, we evaluated the improvements of the AS approach over numerical quadrature schemes in the solution of the EEG-FP. To do so, we utilised a multi-layered spherical head model with anisotropic electrical conductivity field, for which analytical solutions are available [5, 1]. We considered four compartments representing the scalp, skull, cerebrospinal fluid (CSF), and brain, each with outer radius equal to 0.092m, 0.086m, 0.08m, and 0.078m, respectively. The electrical conductivities were considered isotropic for the scalp, CSF, and brain compartments, and set to 0.33 S/m, 1.79 S/m, and 0.33 S/m, respectively. The skull was modelled as an anisotropic layer with radial/tangential conductivities of 0.0093/0.015 S/m. The electric conductivity values were extracted from the relevant literature [26, 27, 28, 29].

We discretised the spherical model using the ISO2Mesh toolbox [30]. Meshes were built to achieve a maximum radius-edge factor of 1.2. We considered a coarser mesh resolution in the brain layer since the volume integral in (2) vanishes in there for sharing the same electrical conductivity as the source neighbourhood. Six discretisations with 39k, 119k, 281k, 440k, 640k, and 938k nodes were obtained. A total of 162 electrodes uniformly placed on the scalp surface were utilised [31]. A visual representation of the models is available in the Supplementary material.

We simulated radially- and tangentially-oriented sources in 100 different locations at a distance $r_{0}$ from the next conductivity jump, with $r_{0}$ taking values in the range $0.125-1.5mm$. For each of them, we computed the relative error (RE), relative difference measure (RDM), and magnification error (MAG), defined as

[TABLE]

with \mbox{\boldmathu}_{n} and \mbox{\boldmathu}_{a} being the numerical and analytical solutions of the EEG-FP, respectively. The RE is a measure of the overall difference between the analytical and numerical solutions of the EEG-FP, whereas the RDM and MAG account for topographical and magnitude differences between them, respectively. Error measures were computed for the solutions obtained with the AS and FS approaches, the latter with $n=2$ and $n=4$.

In addition to standard error measures, we computed the RE between the numerical approximations provided by the AS and FS methods,

[TABLE]

where \mbox{\boldmathu}_{as} and \mbox{\boldmathu}_{fs} are the numerical solutions computed with the AS and FS approaches, respectively. This allowed to evaluate the errors introduced by the utilisation of a numerical quadrature scheme, which are avoided using the AS approach. We calculated the $RE_{s}$ for all mesh discretisations considering the aforementioned sources.

5.3 Illustration in a real model

Finally, we illustrated the advantages of the AS method over the FS approach in a personalised head model. We constructed a realistic head representation based on the Colin 27 high resolution MRI segmentation of the Montreal Neurological Institute [32]. A mesh with 3.8M tetrahedral elements (600k nodes) was created using the ISO2Mesh software as described before. This included the scalp, skull, CSF, grey matter (GM), and white matter (WM) compartments, each of them assumed isotropic and with electrical conductivities equal to 0.33 S/m, 0.004 S/m, 1.79 S/m, 0.45 S/m, and 0.13 S/m, respectively [29]. A total of 162 electrodes were placed on the scalp according to the ABC standard. A visual representation of the model is available in the Supplementary material.

We computed the leadfield matrix for sources on a surface located at the mid-point between the GM/CSF and WM/GM layers, as provided by FreeSurfer [33]. Since the cortical thickness is in the range 1-5mm [34], extra mesh refinement would be needed in thin cortical regions to avoid the FS method to become unstable. We calculated the leadfield matrices with the AS and FS ($n=2$) techniques and obtained the $RE_{s}$. As in the spherical case, this allowed to evaluate the errors introduced in the solution of the EEG-FP due to the use of numerical quadrature.

6 Results

6.1 Local errors

The REs between the element vectors obtained with the analytical and numerical approaches are shown in Fig. 2. Figs. 2a. and b. display the $RE_{e}$ for the surface and volume element source vectors, respectively, as a function of $d/a$ and for integration orders 2, 4, and 6. In addition to providing more accurate solutions, the analytical method allowed to reduce the computational demands on their calculation. In the case of the surface element vector, the analytical method resulted 2, 6, and 12 times faster compared to the numerical technique with $n$ equal to 2, 4, and 6, respectively (in average over 500 iterations). Regarding the volume element vector, the speed-up of the analytical method was 1.4, 4.2, and 11.1 for the same integration orders and number of iterations. It is worth mentioning that the same results were obtained for a given $d/a$ independently of the values chosen for $a$ and $d$, highlighting the importance of such ratio in the error introduced by the FS methodology.

6.2 Spherical head model

In Fig. 3 we show the RE, RDM, and MAG as a function of the distance to the closest conductivity change considering the model with 256k nodes and the AS and FS methods, the latter with $n=2$ and $n=4$. Results are displayed for tangentially-oriented and highly eccentric sources for which the numerical method was shown to introduce noticeable errors. Error measures for radially-oriented sources and other models are presented in the Supplementary material. Regarding the computational cost of the methods, the AS approach resulted, in average, 1.04 and 2.11 times faster than using the FS technique with $n$ equal to 2 and 4, respectively.

Fig. 4 displays the percentiles of the RE obtained for different source locations and model discretisations. In Fig. 4a. we show the percentile curves of the RE considering 100 tangentially-oriented sources located at 0.5mm to the next compartment and all models. As expected, the errors obtained with the AS approach are always lower or equal to those obtained with the FS method regardless of the integration order and head model. Such difference becomes more apparent as the source is closer to the next compartment and the model is less refined. Fig. 4b. presents the 90th percentile of the RE as a function of the distance to the next compartment and the number of first order nodes for both AS and FS ($n=2$ and $n=4$) techniques. It is evident that the FS approach provides solutions with unacceptable REs as the source becomes closer to the closest conductivity jump. Although the increase of the integration order allows to reduce the distance for which the method becomes unstable, this issue cannot be completely avoided by using a numerical quadrature scheme. This is not the case of the AS method, which provides stable solutions for arbitrarily eccentric sources.

The mean RE between the AS and FS approximations is presented in Fig. 5. There, we plot the $RE_{s}$ as a function of the ratio $d/a$, considering all model discretisations and integration orders $n=2$ and $n=4$. We considered $a$ as the mean element side-length belonging to the CSF compartment (i.e. those belonging to the closest compartment with different electrical conductivity). It can be noted that the approximations obtained with the FS method differ from those obtained with the AS approach consistently as a function of $d/a$ regardless of the model discretisation used. In the case of adopting the FS method with $n=2$, errors due to the use of a numerical integration scheme start becoming noticeable (i.e. larger than $1\%$) for $d/a<0.5$. This threshold reduces to $d/a<0.25$ if the integration order is increased to $n=4$.

6.3 Real model

Figs. 6a. and b. present the ratio $d/a$ for the model and sources described in Section 5.3. It can be clearly noted that there are some cortical regions for which the ratio is 0.5 or lower (dark blue in the figures), indicating that the FS ($n=2$) method may introduce unwanted errors compared to the AS approach. This assumption is confirmed in Figs. 6c. and d., which show the $RE_{s}$ between both numerical approximations. It becomes evident that the use of the AS method outperforms the FS technique in these regions, improving up to 7% in the resulting RE.

7 Discussion

We presented analytical expressions for all the matrices and vectors involved in the FE formulation of the subtraction version of the EEG-FP considering dipolar sources. The adoption of these formulas allowed to minimise the errors for a given mesh discretisation, outperforming any subtraction-based approximation available in the literature. We found that the FS approach, considered the most accurate methodology for solving the EEG-FP based on the subtraction formulation [13], introduces errors depending on the mesh refinement and source location through the ratio $d/a$. These inaccuracies were minimised with the AS technique, which prevents errors growing disproportionately no matter the mesh density and source position. Moreover, these advantages come at no extra computational cost, which resulted slightly lower than that corresponding to the FS method with $n=2$. This makes the AS technique the optimal option for solving the EEG-FP for any given model discretisation with regard to accuracy.

The experiments presented provide new insights into the nature of the local errors due to the use of numerical quadrature formulas, as well as their impact in the overall solution of the EEG-FP. In the case of the surface element vector, $d$ (i.e. the distance between the source and the scalp surface) is always larger than 14mm. Then, from Fig. 2a., errors in \mbox{\boldmathb}^{s} introduced by a numerical scheme would become negligible for surface triangles with average side-length smaller than 14mm, which is always the case. This leaves the approximation of \mbox{\boldmathb}^{v} as the main responsible for the increased errors in the EEG-FP by means of the FS approach compared to those obtained with the AS method. In this case, errors will appear in the integration over elements with electrical conductivity different from that of the source compartment. We found that the use of numerical quadrature introduces errors depending on the ratio $d/a$ and the integration order $n$ (Figs. 2 and 5). These errors became noticeable for $d/a<0.5\ (0.25)$ in case of using $n=2\ (4)$, reaching values of $RE=1$ for $d/a\approx 0.05$. Such errors are avoided with the AS technique, which provides stable results regardless of $d/a$.

Based on a real head model, we showed that the use of the AS method has clear benefits compared to the use of the FS approach, specially for sources located in the somatosensory and visual cortices, which are known for having an average thickness of approximately 1mm [34]. Considering sources located in the midpoint between the GM/CSF and WM/GM surface layers, whose distance to the CSF and WM compartments is in the range 0.5-2.25mm [34], a tetrahedral mesh with side-lengths smaller than 1mm would be needed to avoid inaccuracies with the FS methodology (for $n=2$). Such refinement would need to take place in both CSF and WM regions, resulting in a highly-refined model and, consequently, in very high computational requirements. On the contrary, the AS method would not need such refinement to lead to even more accurate results, reducing both errors and computational demands. In the model used in Section 5.3, this resulted in the FS approach introducing an additional 7% of RE in selected brain regions compared to the AS method, highlighting the value of the latter.

Detailed head models are being increasingly considered for the analysis of both invasive and non-invasive EEG recordings. Recent results have demonstrated the value of incorporating the vessel network [35] and even the dura [36] into the models used to localise sources of electrical activity. In this direction, evidence point to the incorporation of inhomogeneity of the electrical conductivity profile of the brain cortex, which presents significant laminar dependencies with radial/tangential anisotropies [37]. The addition of such level of detail would require a numerical method stable for extremely-eccentric sources located very close to conductivity jumps, being the AS approach the most appropriate (if not the only) option.

Within all the methodologies available for simulating electromagnetic fields generated by brain current sources, the boundary element method is by far the most utilised. One of the reasons for this phenomenon is based on the great levels of accuracy reported for isotropic layered models, which are partly supported by the existence of analytical formulas for the potential integrals involved [38, 39]. Such expressions were shown fundamental for speeding up simulations and avoiding errors due to the utilisation of a numerical quadrature scheme. We believe that the AS approach will have a similar effect with the FEM, extending its adoption and allowing optimal results in anisotropic and highly-detailed head representations.

The MATLAB implementation of the method is publicly available through the FEMEG toolbox (https://femeg.github.io).

8 Conclusions

We presented analytical expressions for all the integrals required by the FE formulation of the EEG-FP based on the subtraction approach. Such formulas were shown to increase the accuracy of the solution of the EEG-FP, preventing errors growing disproportionately independently of the mesh density and source position (unlike any other existing method). We demonstrated that the computational cost is similar to that obtained with the lowest order numerical scheme, making the AS method a competitive option in the field. We believe this approach will be considered as a gold standard technique for solving the EEG-FP considering fully-realistic head models.

Appendix 1: Proof of Proposition 1

We use [24, A3.31] to rewrite the surface divergence,

[TABLE]

The final result is obtained by noting that

[TABLE]

Appendix 2: Proof of Proposition 2

For any two vector fields $a$ and $b$ in $\mathbb{R}^{3}$, it is easily verified that the following relation holds

[TABLE]

where \nabla\mbox{\boldmatha} is the gradient of the vector field with elements \nabla\mbox{\boldmatha}_{ij}=\partial a_{i}/\partial x_{j}. Let \mbox{\boldmatha}=\mbox{\boldmathq} and \mbox{\boldmathb}=\mbox{\boldmathR}R^{-3}. Direct computation of the partial derivatives of \mbox{\boldmathR}R^{-3} allows to check that its gradient is symmetric, i.e. \nabla(\mbox{\boldmathR}R^{-3})=\nabla(\mbox{\boldmathR}R^{-3})^{T}. Since $q$ is constant, the application of (32) yields

[TABLE]

Appendix 3: Derivation of (28)

We operate as described in Section 4.3, i.e. split the integral over $T$ as the sum of the integrals over $T_{\epsilon}$ and $T-T_{\epsilon}$, after which we take the limit $\epsilon\rightarrow 0$. Utilising (15) we get

[TABLE]

Let \alpha(\mbox{\boldmath\rho}) be the portion of the circular arc $T_{\epsilon}$ lying within $T$. This parameter is equal to [math] if $\rho$ is outside $T$, $2\pi$ if $\rho$ is inside $T$, $\pi$ if $\rho$ is on $\partial T$, or $0<\alpha<\pi$ if it is on a node. Utilising polar coordinates, the first two integrals yield

[TABLE]

and

[TABLE]

The computation of the third integral in (33) is performed in $\mathcal{L}_{a}$. Recalling the parameterisation \mbox{\boldmaths}(s)=\mbox{\boldmatht}_{i}^{0}+s\hat{\mbox{\boldmaths}}_{i}, valid for \mbox{\boldmaths}\in\partial T, and noting that \hat{\mbox{\boldmathm}}_{i}\cdot\mbox{\boldmatht}=t_{i}^{0} for \mbox{\boldmatht}\in\partial T, we get

[TABLE]

The final step consists of replacing these expressions into (33) and take the limit $\epsilon\rightarrow 0$. This leads to

[TABLE]

Noting that \alpha(\mbox{\boldmath\rho})=\sum_{i=1}^{3}\beta_{i}(w_{0}=0) [23, eq. (17)], and employing the identity [18]

[TABLE]

we get the desired result.

Appendix 4: Derivation of (30)

We multiply both sides of (25) by $\{u_{a},v_{a}\}^{T}$ and integrate over $T$,

[TABLE]

After utilising (18), the first term yields

[TABLE]

The integrals in (35) are solved in $\mathcal{L}_{a}$. Using (12) we get

[TABLE]

where we used

[TABLE]

The second integral in (34) is solved using (16), resulting in

[TABLE]

Similarly, we apply (17) to solve the third integral in (34), obtaining

[TABLE]

Finally, we replace (27) and (36)–(38) into (34), leading to the expected result.