Handling Anisotropic Conductivities in the EEG Forward Problem with a Symmetric Formulation

TL;DR

This paper introduces a modified symmetric surface integral formulation for the EEG forward problem that effectively models anisotropic conductivities in the skull and white matter, improving accuracy in source analysis.

Contribution

It presents a novel surface integral equation approach that incorporates anisotropic conductivity profiles using an anisotropic fundamental solution.

Findings

The proposed method accurately models anisotropic conductivities.

Numerical results demonstrate improved forward problem solutions.

Impact shown on both canonical and real EEG scenarios.

Abstract

The electroencephalography (EEG) forward problem, the computation of the electric potential generated by a known electric current source configuration in the brain, is a key step of EEG source analysis. In this problem, it is often desired to model the anisotropic conductivity profiles of the skull and of the white matter. These profiles, however, cannot be handled by standard surface integral formulations and the use of volume finite elements is required. Leveraging on the representation theorem using an anisotropic fundamental solution, this paper proposes a modified symmetric formulation for solving the EEG forward problem by a surface integral equation which can take into account anisotropic conductivity profiles. A set of numerical results is presented to corroborate theoretical treatments and to show the impact of the proposed approach on both canonical and real case scenarios.

Tables

Table 1. Table 1: Timing results to obtain an error of 1 % percent 1 1\% with respect to the analytical anisotropic solution. ”Time Z” is the time to build the system matrix, ”Time sol” is the time to solve the system (direct inversion), ”N elements”is the number of basis functions, sparsity is the ratio of the number of non-zero of the matrix with respect to the size of the matrix.

Time Z	Time sol.	N elements	sparsity	Relative Error
$545$ s	$12.7$ s	$19,392$	$0.1$	$0.013$

Equations

$${\gamma }_{0i}^{\pm }g=g_{|{\Gamma }_i^{\pm }}$$

$${\gamma }_{1_i}^{\pm }g=\vec{n}\cdot \overline{\stackrel{ˉ}{\sigma }}\nabla g_{|{\Gamma }_i^{\pm }}.$$

$$\nabla \cdot \overline{\stackrel{ˉ}{\sigma }}\nabla V=f$$

$${\left[V\right]}_k=0$$

$${\left[\vec{n}\cdot \overline{\stackrel{ˉ}{\sigma }}\nabla V\right]}_k=0.$$

$$G_j(r,r^{\prime })=\frac{1}{4\pi }\frac{1}{\sqrt{det({\overline{\stackrel{ˉ}{\sigma }}}_j)}\sqrt{{\overline{\stackrel{ˉ}{\sigma }}}_j^{-1}(r-r^{\prime })\cdot (r-r^{\prime })}}$$

$$V_i(r)={\int }_{{\Omega }_i}{G}_i(r,r^{\prime })f(r^{\prime })d{r}^{\prime }+{\int }_{{\partial }_{{\Omega }_i}}{\overline{\stackrel{ˉ}{\sigma }}}_j\overrightarrow{n^{\prime }}\cdot \nabla G_i(r,r^{\prime })V(r^{\prime })d{r}^{\prime }-{\int }_{{\partial }_{{\Omega }_i}}{\overline{\stackrel{ˉ}{\sigma }}}_j\overrightarrow{n^{\prime }}\cdot \nabla V(r^{\prime })G_i(r,r^{\prime })d{r}^{\prime }.$$

$${\mathcal{S}}_i\nu (r)={\int }_{{\partial }_{{\Omega }_i}}{G}_i(r,r^{\prime })\nu (r^{\prime })d{r}^{\prime }{\mathcal{S}}_i:H^{-1/2}({\partial }_{{\Omega }_i})\to H^1(\Omega )$$

$${\mathcal{D}}_i\mu (r)={\int }_{{\partial }_{{\Omega }_i}}\mu (r^{\prime })\overrightarrow{n^{\prime }}.(\overline{\stackrel{ˉ}{{\sigma }_i}}\nabla G_i(r,r^{\prime }))d{r}^{\prime }{\mathcal{D}}_i:H^{1/2}({\partial }_{{\Omega }_i})\to H^1(\Omega ).$$

$$V_i(r)={\int }_{{\Omega }_i}{G}_i(r,r^{\prime })f(r^{\prime })d{r}^{\prime }+{\mathcal{D}}_{i}V(r)-{\mathcal{S}}_i\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i\nabla V(r^{\prime }).$$

$${\gamma }_{0i}^{\pm }{\mathcal{S}}_j=\sum _{k\in I_j}{S}_{ik,{\Omega }_j}$$

$${\gamma }_{1i}^{\pm }{\mathcal{S}}_j=\mp \frac{1}{2}+\sum _{k\in I_j}{D}_{ij,{\Omega }_j}^{\ast }$$

$${\gamma }_{0i}^{\pm }{\mathcal{D}}_j=\mp \frac{1}{2}+\sum _{k\in I_j}{D}_{ij,{\Omega }_j}$$

$${\gamma }_{1i}^{\pm }{\mathcal{D}}_j=\sum _{k\in I_j}{N}_{ij,{\Omega }_j}$$

$$S_{ik,{\Omega }_j}\nu (r)={\int }_{{\Gamma }_k}{G}_j(r,r^{\prime })\nu (r^{\prime })d{r}^{\prime },S_{ij}:H^{-1/2}({\Gamma }_k)\to H^{1/2}({\Gamma }_i),$$

$$D_{ik,{\Omega }_j}^{\ast }\nu (r)={\int }_{{\Gamma }_k}\nu (r^{\prime })\vec{n}\cdot ({\overline{\stackrel{ˉ}{\sigma }}}_j\nabla G_j(r,r^{\prime }))d{r}^{\prime },D_{ij}^{\ast }:H^{-1/2}({\Gamma }_k)\to H^{-1/2}({\Gamma }_i),$$

$$D_{ik,{\Omega }_j}\mu (r)={\int }_{{\Gamma }_k}\mu (r^{\prime })\overrightarrow{n^{\prime }}\cdot ({\overline{\stackrel{ˉ}{\sigma }}}_j\nabla G_j(r,r^{\prime }))d{r}^{\prime },D_{ij}:H^{1/2}({\Gamma }_k)\to H^{1/2}({\Gamma }_i),\text{ and }$$

$$N_{ik,{\Omega }_j}\mu (r)=\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i\nabla {\int }_{{\Gamma }_k}\mu (r^{\prime })\overrightarrow{n^{\prime }}\cdot ({\overline{\stackrel{ˉ}{\sigma }}}_j\nabla G_j(r,r^{\prime }))d{r}^{\prime },N_{ij}:H^{1/2}({\Gamma }_k)\to H^{-1/2}({\Gamma }_i).$$

$$V_j(r)={\int }_{{\Omega }_j}{G}_j(r,r^{\prime })f(r^{\prime })d{r}^{\prime }+{\mathcal{D}}_{j}V(r)-{\mathcal{S}}_j\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_j\nabla V(r).$$

$${\gamma }_{0_k}^{+}{V}_j={\gamma }_{0_k}^{+}{v}_{di{p}_j}+\sum _{l\in I_j}{D}_{kl,{\Omega }_j}{\mathcal{V}}_l+\frac{1}{2}{\mathcal{V}}_k+\sum _{l\in I_j}{S}_{kl,{\Omega }_j}\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_j\nabla {\mathcal{V}}_l.$$

$${\gamma }_{0_k}^{-}{V}_i={\gamma }_{0_k}^{-}{v}_{di{p}_i}+\sum _{l\in I_i}{D}_{kl,{\Omega }_i}{\mathcal{V}}_l+\frac{1}{2}{\mathcal{V}}_k-\sum _{l\in I_i}{S}_{kl,{\Omega }_i}\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i\nabla {\mathcal{V}}_l.$$

$$\begin{array}{rl}\sum _{l\in I_i}{D}_{kl,{\Omega }_i}{\mathcal{V}}_l-\sum _{l\in {\omega }_i}{D}_{kl,{\Omega }_j}{\mathcal{V}}_l& +\sum _{l\in I_i}{S}_{kl,{\Omega }_i}\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i\nabla {\mathcal{V}}_l\\ & -\sum _{l\in I_j}{S}_{kl,{\Omega }_j}\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_j\nabla {\mathcal{V}}_l=-\left({\gamma }_{0_k}^{-}{v}_{di{p}_i}-{\gamma }_{0_k}^{+}{v}_{di{p}_j}\right)\end{array}$$

$${\gamma }_{1_k}^{-}{V}_i={\gamma }_{1_k}^{-}{v}_{di{p}_i}+\sum _{l\in I_i}{N}_{kl,{\Omega }_i}{\mathcal{V}}_l+\frac{1}{2}\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i{\mathcal{V}}_k-\sum _{l\in I_i}{D}_{kl,{\Omega }_i}^{\ast }\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i\nabla {\mathcal{V}}_l$$

$${\gamma }_{1_k}^{+}{V}_j={\gamma }_{1_k}^{+}{v}_{di{p}_j}+\sum _{l\in I_j}{N}_{kl,{\Omega }_j}{\mathcal{V}}_l+\frac{1}{2}\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_j{\mathcal{V}}_k-\sum _{l\in I_j}{D}_{kl,{\Omega }_j}^{\ast }\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_j\nabla {\mathcal{V}}_l.$$

$$\begin{array}{rl}\sum _{l\in I_i}{N}_{kl,{\Omega }_i}{\mathcal{V}}_l-\sum _{l\in I_j}{N}_{kl,{\Omega }_j}{\mathcal{V}}_l& +\sum _{l\in I_i}{D}_{kl,{\Omega }_i}^{\ast }\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_i\nabla {\mathcal{V}}_l\\ & -\sum _{l\in I_j}{D}_{kl,{\Omega }_j}^{\ast }\vec{n}\cdot {\overline{\stackrel{ˉ}{\sigma }}}_j\nabla {\mathcal{V}}_l=-\left({\gamma }_{1_k}^{-}{v}_{di{p}_i}-{\gamma }_{1_k}^{+}{v}_{di{p}_j}\right).\end{array}$$

$$\mathbf{Za}=\mathbf{b}$$

$${\left[\mathbf{Z}\right]}_{kl}={\int }_{{\mu }_k^{\ast }}{f_t}_k(r^{\prime })Z{f_e}_l(r^{\prime })d{r}^{\prime }$$

$${\left[{\mathbf{S}}_{\mathbf{ij},{\mathbf{\Omega }}_{\mathbf{k}}}\right]}_{kl}={\int }_{{\mu }_{P_{0k}}}{P}_{0k}{S}_{ij,{\Omega }_k}{P}_{0l}(r)dr$$

$${\left[{\mathbf{D}}_{\mathbf{ij},{\mathbf{\Omega }}_{\mathbf{k}}}\right]}_{kl}={\int }_{{\mu }_{P_{0k}}}{P}_{0k}{D}_{ij,{\Omega }_k}{P}_{1l}(r)dr$$

$${\left[{\mathbf{D}}_{\mathbf{ij},{\mathbf{\Omega }}_{\mathbf{k}}}^{\ast }\right]}_{kl}={\int }_{{\mu }_{P_{1k}}}{P}_{1k}{D}_{ij,{\Omega }_k}^{\ast }{P}_{0l}(r)dr$$

$${\left[{\mathbf{N}}_{\mathbf{ij},{\mathbf{\Omega }}_{\mathbf{k}}}\right]}_{kl}={\int }_{{\mu }_{P_{1k}}}{P}_{1k}{N}_{ij,{\Omega }_k}{P}_{1l}(r)dr.$$

$$\mathbf{Z}=\left[\begin{array}{cc}\mathbf{S}& \mathbf{D}\\ {\mathbf{D}}^{\ast }& \mathbf{N}\end{array}\right],$$

$$\mathbf{A}=\left(\begin{array}{cccccc}{\mathbf{A}}_{11,{\mathbf{\Omega }}_1}-{\mathbf{A}}_{11,{\mathbf{\Omega }}_2}& {\mathbf{A}}_{12,{\mathbf{\Omega }}_1}& -{\mathbf{A}}_{13,{\mathbf{\Omega }}_2}& -{\mathbf{A}}_{14,{\mathbf{\Omega }}_2}& 0& 0\\ {\mathbf{A}}_{21,{\mathbf{\Omega }}_1}& {\mathbf{A}}_{22,{\mathbf{\Omega }}_1}-{\mathbf{A}}_{22,{\mathbf{\Omega }}_3}& -{\mathbf{A}}_{23,{\mathbf{\Omega }}_3}& 0& -{\mathbf{A}}_{25,{\mathbf{\Omega }}_3}& 0\\ {\mathbf{A}}_{31,{\mathbf{\Omega }}_2}& -{\mathbf{A}}_{32,{\mathbf{\Omega }}_3}& {\mathbf{A}}_{33,{\mathbf{\Omega }}_2}-{\mathbf{A}}_{33,{\mathbf{\Omega }}_3}& {\mathbf{A}}_{34,{\mathbf{\Omega }}_2}& -{\mathbf{A}}_{35,{\mathbf{\Omega }}_3}& 0\\ {\mathbf{A}}_{41,{\mathbf{\Omega }}_2}& 0& {\mathbf{A}}_{43,{\mathbf{\Omega }}_2}& {\mathbf{A}}_{44,{\mathbf{\Omega }}_2}-{\mathbf{A}}_{44,{\mathbf{\Omega }}_4}& -{\mathbf{A}}_{45,{\mathbf{\Omega }}_4}& -{\mathbf{A}}_{46,{\mathbf{\Omega }}_4}\\ 0& {\mathbf{A}}_{52,{\mathbf{\Omega }}_3}& {\mathbf{A}}_{53,{\mathbf{\Omega }}_3}& -{\mathbf{A}}_{54,{\mathbf{\Omega }}_4}& {\mathbf{A}}_{55,{\mathbf{\Omega }}_3}-{\mathbf{A}}_{55,{\mathbf{\Omega }}_4}& -{\mathbf{A}}_{56,{\mathbf{\Omega }}_4}\\ 0& 0& 0& {\mathbf{A}}_{64,{\mathbf{\Omega }}_4}& {\mathbf{A}}_{65,{\mathbf{\Omega }}_4}& {\mathbf{A}}_{66,{\mathbf{\Omega }}_4}\end{array}\right),$$

$${\left[{\mathbf{b}}_{{\gamma }_{\mathbf{p}}\mathbf{i}}\right]}_k={\int }_{{\mu }_k}{f_t}_k(r^{\prime })b_{{\gamma }_{p}i}(r^{\prime })d{r}^{\prime },$$

$$\{\begin{array}{rl}& {\left[{\mathbf{a}}_{0j}\right]}_l={\beta }_{j,l}\\ & {\left[{\mathbf{a}}_{1j}\right]}_l={\alpha }_{j,l}\end{array}$$

Full text

††journal: Physics in Medicine and Biology

1 Introduction

Epileptic source localization from high resolution electroencephalographies (EEGs) is a fundamental step in the pre-surgical evaluation of focal epileptic patients that are refractory to pharmacological treatments and for whom surgical resection of the epileptic focus is considered [1, 2, 3, 4]. In this brain imaging technique, starting from the electrical potential measured on the scalp, the brain current sources responsible for focal epilepsy are localized by solving an inverse source problem [5, 6, 7]. Solving this inverse problem requires the multiple solutions of an EEG forward problem that provides, from known brain electrical current sources, the surface potential measured at the electrodes’ locations [5, 8]. The solution of the forward problem must be computed with the highest possible precision [9, 10, 11] to increase the accuracy of the source localisation process.

In solving the EEG forward problem, spherical head models have been historically used since analytic solutions are available for them [12, 13]. However, modern techniques rely on the use of realistic head models that require a numerical solution, but for which the accuracy of the forward EEG solution is largely improved [14, 15, 16, 17, 18]. Several methods can be used for numerically solving the EEG forward problem, including Finite Difference (FDM) [19, 20] , Finite Element (FEM) [21] and Boundary Element methods (BEM) [22, 23, 24, 25]. Boundary Element methods have been quite popular given that they require only surfacic discretization of the brain layers when compared to the other two techniques that rely on volume discretizations. These methods however are not panacea given that, in their standard incarnations, they cannot handle anisotropic conductivity profiles. Indeed, correct modeling of anisotropic conductivity profiles is quite important given the influence of white matter and skull anisotropic conductivities on source localization [26, 27, 28, 29, 30, 31, 32]. The reader should notice that anisotropic conductivity profiles can be naturally treated with FEM approaches. Even if the computational load can be reduced by resorting to transfer matrices, as in [33] or [34], a purely surfacic BEM method (with anisotropic modeling capabilities) could be desirable. Indeed, in solving the EEG forward problem, the integral methods, which are based on an operator having an asymptotically smooth kernel, can easily be used with a fast multipole method (FMM) [35]. This decreases their complexity to a quasi-linear one ($O(N\log(N))$, where $N$ is the number of unknowns). Moreover, exploiting the reciprocity theorem as in [19, 36], further reduces the computational burden when solving the EEG inverse problem. Some steps in building an anisotropic BEM formulation have been presented in [37] where virtual sources are placed in the fibers to account for their anisotropic properties. Moreover, [38] proposed an interesting coordinate transform to handle constant anisotropy in a single conducting body. No details, however, were provided for the multi-compartments problem and the associated numerical solutions. Preliminary and partial results on an indirect anisotropic formulation for the EEG forward problem had been presented in the conference contributions [39, 40].

This work proposes a new approach to the problem that falls under the family of symmetric formulations as originally introduced in [24] and results in a surfacic discretized BEM formulation capable of handling constant piecewise homogeneous conductivity profiles including anisotropies for both nested and non-nested compartments. This is obtained by leveraging on the representation theorem for the Poisson equation, using an anisotropic fundamental solution. Instead of building a harmonic solution as is commonly done in building EEG integral formulations, our approach directly computes the unknown potential. This allows to enforce boundary conditions also across different anisotropic media. This approach also naturally permits to take into account non-nested domains as in [41], and it therefore renders possible further subdivisions of the compartments according to the variation of their conductivity tensors.

The article is organized as follows: Section 2 describes background material and sets the notations. Section 3 presents the new anisotropic symmetric integral equation. Section 4 focuses on the discretization of the BEM solution and its implementation. Section 5 complements all the theoretical developments with numerical results that show the effectiveness of the newly proposed method in both canonical and real case scenarios.

2 Background and Notation

Consider a domain $\Omega$ divided into $N$ non overlapping subdomains $\Omega_{i}$ with Lipschitz boundary $\partial_{\Omega_{i}}$ such that $\Omega=\bigcup_{j}^{N}\overline{\Omega}_{i}$. The external domain, $\mathbb{R}^{3}\backslash\Omega$, is denoted with $\Omega_{N+1}$. We further define $\Gamma_{ij}=\overline{\Omega}_{i}\cap\overline{\Omega}_{j}$ the interface beween $\Omega_{i}$ and $\Omega_{j}$ as well as $\omega_{i}$ the set of neighboring domains of $\Omega_{i}$ so that $\omega_{i}=\{j\leq N+1|\overline{\Omega}_{i}\cap\overline{\Omega}_{j}\neq\emptyset\}$. Then $\partial_{\Omega_{i}}$ can be defined as $\partial_{\Omega_{i}}=\bigcup_{k\in I_{i}}\Gamma_{k}$ with $I_{i}=\{k|\Gamma_{k}=\Gamma_{ij},j\in\omega_{i}\}$, the set of interfaces composing $\partial_{\Omega_{i}}$. The normal $\vec{n}_{ij}$ to the interface $\Gamma_{ij}$ is oriented from $\Omega_{i}$ to $\Omega_{j}$ with $i<j$. Figure 1 shows a general decomposition of $\Omega$ into subdomains. Standard EEG BEM formulations often make use of nested domains [24], so that when a three layers geometry is chosen, the compartments represents the brain, the skull and the scalp. In this particular case, $\partial_{\Omega_{i}}=\Gamma_{i,i-1}\cup\Gamma_{i,i+1}$. We define the traces of a function $g$ on a boundary $\Gamma_{i}$ and of its conormal derivative as [42]

[TABLE]

where $\bar{\bar{\sigma}}$ is the conductivity tensor. Moreover $[\cdot]_{k}$ will denote the jump of a function across the surface $\Gamma_{k}$: $[g]_{k}=\gamma_{0k}^{-}g_{k}-\gamma_{0k}^{+}g_{k}$ and $\left[\vec{n}\cdot\bar{\bar{\sigma}}\nabla g\right]_{k}=\gamma_{1k}^{-}g_{k}-\gamma_{1k}^{+}g_{k}$.

The EEG forward problem amounts at computing the electric potential $V$ knowing the brain electric sources $f=\nabla\cdot j$ when the current $j$ propagates in a medium of conductivity $\bar{\bar{\sigma}}$ which is a real symmetric and positive definite matrix. Under standard quasi-static assumptions this calls for the solution of the Poisson’s equation [43]

[TABLE]

with boundary conditions at each interface $\Gamma_{k}$

[TABLE]

The conditions above enforce the continuity of the potential and of its derivative between different compartments of the head. The conductivity $\bar{\bar{\sigma}}$ is assumed to be piecewise homogeneous and potentially anisotropic, $\bar{\bar{\sigma}}_{i}$ will denote the conductivity of the domain $\Omega_{i}$. The source term $f=\nabla\cdot j$ is usually a linear combination of dipole sources $f_{j}$ such that $f_{j}=q_{j}\cdot\nabla\delta_{r_{0j}}$ with $q_{j}$ the dipole moment and $r_{0j}$ the dipole position [44], [45].

3 A Surface Integral Formulation for Anisotropic Conductivity Profiles

The fundamental solution (Green’s function) of (2) in an unbounded medium reads [46]

[TABLE]

where $\bar{\bar{\sigma}}_{j}$ denotes the homogeneous conductivity tensor of the domain $\Omega_{j}$. Decomposing the potential $V$ into different contributions in each domains $\Omega_{i}$, $V=\sum_{i=1}^{N}V_{i}$ such that $V_{i}(r)=V(r)$ if $r\in\Omega_{i}$ and 0 elsewhere, the representation theorem in $\Omega_{i}$ gives [42]

[TABLE]

Define the following integral operators

[TABLE]

and

[TABLE]

Here and in the following $H^{s}$ will denote the Sobolev space of order $s$. Then, since $\bar{\bar{\sigma_{i}}}$ is symmetric, the representation theorem (5) reads

[TABLE]

Note that, with our notations and normal orientation, [42]

$\displaystyle\gamma_{0i}^{\pm}\mathcal{S}_{j}=\sum_{k\in I_{j}}S_{ik,\Omega_{j}}$

(9a)

$\displaystyle\gamma_{1i}^{\pm}\mathcal{S}_{j}=\mp\frac{1}{2}+\sum_{k\in I_{j}}D^{*}_{ij,\Omega_{j}}$

(9b)

and

$\displaystyle\gamma_{0i}^{\pm}\mathcal{D}_{j}=\mp\frac{1}{2}+\sum_{k\in I_{j}}D_{ij,\Omega_{j}}$

(10a)

$\displaystyle\gamma_{1i}^{\pm}\mathcal{D}_{j}=\sum_{k\in I_{j}}N_{ij,\Omega_{j}}$

(10b)

with,

[TABLE]

Now, let us consider $\Gamma_{k}=\overline{\Omega_{i}}\cap\overline{\Omega_{j}}\neq\emptyset$ with $i<j$ to express the interior problems and combine the obtained equations to enforce the boundary conditions at the interface $\Gamma_{k}$ (see figures 2(a), 2(b), and 2(c)). The first step is to express also the interior problem (8) in $\Omega_{j}$

[TABLE]

Then, decomposing $V_{j_{|\partial\Omega_{j}}}$ and $V_{i_{|\partial\Omega_{i}}}$ into their contributions on the different interfaces $\Gamma_{l}$ with $l\in I_{j}$ or $I_{i}$ respectively, such that $V_{j_{|\Gamma_{l}}}=\mathcal{V}_{l}$ for $l\in I_{j}$ or $V_{i_{|\Gamma_{l}}}=\mathcal{V}_{l}$ for $l\in I_{i}$ respectively, and taking the trace of (12) and of (8) we obtain

[TABLE]

[TABLE]

where to simplify the notations we introduced $v_{dip_{j}}(r)=\int_{\Omega_{j}}G_{j}(r,r^{\prime})f(r^{\prime})dr^{\prime}$ the potential due to the source in $\Omega_{j}$. Finally, enforcing the boundary condition (3a) by subtracting (13) to (14) we obtain

[TABLE]

The same reasoning using the trace of the conormal derivative gives the two following interior problems in $\Omega_{j}$ and $\Omega_{i}$

[TABLE]

and

[TABLE]

Given that at the interface $\Gamma_{k}$, $\frac{1}{2}\vec{n}\cdot\bar{\bar{\sigma}}_{i}\mathcal{V}_{k}=\frac{1}{2}\vec{n}\cdot\bar{\bar{\sigma}}_{j}\mathcal{V}_{k}$, the equations (16) and (17) can be combined using the second boundary condition (3b) to get a second integral equation

[TABLE]

The equations (15) and (18) correspond to the symmetric formulation presented in [24]. However, the introduced formulation differs from it in its essence as the kernel of the involved integral operators is not the same. Moreover, the presenting reasoning generalizes the use of this formulation to fully inhomogeneous and anisotropic domains. As in the case of the standard symmetric formulation, at the external interface(s), only the boundary condition (3b) is necessary to be enforced, this means that for this(ese) interface(s), only (18) will be enforced. Indeed, (3b) at the outermost interface directly tells that the conormal derivative of the potential is zero here.

4 Discretization and Implementation of the Modified Symmetric Equation

To obtain a linear system to be solved, following the usual BEM strategy [42], the integral equations must be tested with a suitably chosen set of basis functions. Tessellating the geometry into $N_{t}$ triangular cells $t_{i}$ and $N_{v}$ vertices $v_{i}$, this gives rise to the linear system

[TABLE]

where $\mathbf{Z}$ contains the discrete version of the integral equation. Each entry of $\mathbf{Z}$ is given by

[TABLE]

where $Z$ is the operator of the integral equation under consideration, $\{f_{t_{k}}\}$ is the set of testing functions, and $\{f_{e_{k}}\}$ is the set of expansion functions whose supports are given by $\{\mu_{e_{k}}\}$. In our case, the unknowns $\mathcal{V}_{i}$ and $\vec{n}\cdot\bar{\bar{\sigma_{j}}}\nabla\mathcal{V}_{i}$ can be respectively discretized using, as expansion functions, piecewise linear functions $P_{1k}\in H^{1/2}$ and piecewise constant functions $P_{0k}\in H^{-1/2}$. Such that,

$\mathcal{V}_{j}\approx\sum_{k}\alpha_{j,k}P_{1k}$ and $\bar{\bar{\sigma_{i}}}\nabla\mathcal{V}_{j}\approx\sum_{k}\beta_{j,k}P_{0k}$,

where $P_{0i}(r)=1$ if $r\in t_{i}$ and [math] elsewhere and $P_{1i}(r)=1$ if $r=v_{i}$ and linearly decreases to [math] on their support $\mu_{P_{1i}}=\{t_{k}|v_{i}\textnormal{ is a vertex of }t_{k}\}$. The support of $P_{0i}$ is given by $\mu_{P_{0i}}=t_{i}$.

Expansion and testing functions should comply with the operatorial mappings of the involved integral operators. The range of $S_{ij}$ and $D_{ij}$ is $H^{1/2}$ for which the dual space is $H^{-1/2}$. As $P_{0}\in H^{-1/2}$ these functions can be used for testing equation (15). The operators involved in (18) are $D_{ij}^{*}$ and $N_{ij}$ whose range is $H^{-1/2}$. This means that pyramidal functions $P_{1}\in H^{1/2}$ are suitable testing functions. This discretization corresponds to the one presented in [24]. Then the entry of the matrix $\mathbf{Z}$ will be given by the discrete version of the operators $S_{ij,\Omega_{k}},D_{ij,\Omega_{k}},D^{*}_{ij,\Omega_{k}}$ and $N_{ij,\Omega_{k}}$ with

[TABLE]

Note here that only surface basis functions are required. As such, only a surface mesh is necessary to implement the presented formulation. The separation of the considered geometry into different homogeneous anisotropic or isotropic compartments is however a necessary pre-processing step.

The proposed formulation can be interesting to model the skull anisotropy where in this layer, different conductivity tensors would correspond to different domains. Following this, as well as to fix the ideas and to simplify the implementation, we provide below the explicit form of the matrices in the case of the geometry depicted figure 3. In this simplified case, the system matrix reads

[TABLE]

with

[TABLE]

where $\mathbf{A}$ stands for the operator $\mathbf{S}$, $\mathbf{D}$, $\mathbf{D^{*}}$ or $\mathbf{N}$. Note that the last line of $\mathbf{A}$ must not be computed when $\mathbf{A}=\mathbf{S}$ or $\mathbf{D}$, nor its last column when $\mathbf{A}=\mathbf{S}$ or $\mathbf{D^{*}}$. Indeed, equation (15) is not necessary at the outermost interface and (3b) at the outermost interface directly gives $\vec{n}\cdot\bar{\bar{\sigma_{4}}}\nabla V_{4_{6}}=0$.

The right hand side $\mathbf{b}=\left[\mathbf{b_{\gamma_{0}}}\,\mathbf{b_{\gamma_{1}}}\right]^{T}$ is defined using $\mathbf{b_{\gamma_{p}}}=\left[\mathbf{b_{\gamma_{p}1}}\mathbf{b_{\gamma_{p}2}},...,\mathbf{b_{\gamma_{p}N}}\right]^{T}$ with $p=0$ or $1$, $N=5$ if $p=0$ or $N=6$ if $p=1$, and

[TABLE]

where $b_{\gamma_{p}i}(r)=-\left(\gamma_{p_{i}}^{-}v_{dip_{k}}(r)-\gamma_{p_{i}}^{+}v_{dip_{l}}(r)\right)$, and ${f_{t}}_{k}=P_{1k}$ if $p=0$ and ${f_{t}}_{k}=P_{0k}$ if $p=1$. In the same fashion, the vector $\mathbf{a}$ that contains the coefficients in the unknown expansion reads $\mathbf{a}=\left[\mathbf{a_{0}}\,\mathbf{a_{1}}\right]^{T}$ with $\mathbf{a_{p}}=\left[\mathbf{a_{p1}},\mathbf{a_{p2}},...,\mathbf{a_{pN}}\right]^{T}$ with $p=0$ or $1$, $N=5$ if $p=0$ or $N=6$ if $p=1$, and

[TABLE]

4.1 Computation of the innermost integral

In practice, computing the elements $[\mathbf{A_{i}}]_{kl}$ requires two integrations. The inner integral requires to compute either the integral of $G_{j}$ or of its derivative times the expansion function $f_{el}$. An analytical solution to obtain this integration for the usual Green’s function $G(r,r^{\prime})=\frac{1}{4\pi}\frac{1}{||r-r^{\prime}||}$ has been proposed in [47]. In (4), the change of variable $R=\sqrt{\bar{\bar{\sigma}}_{i}^{-1}}r$ (${\bar{\bar{\sigma}}}_{i}$ is a symmetric positive definite tensor) transforms $G_{i}(r,r^{\prime})$ into $G_{i}(R,R^{\prime})=\frac{1}{4\pi}\frac{1}{||R-R^{\prime}||}$ for which we can apply the analytical integration formulas in [47]. The outer integrals (arising from the testing of the operators) are performed numerically using Gaussian integration rules.

4.2 Computation of the hypersingular operator

Using the following representation [42]

[TABLE]

$\mathbf{N_{ij,\Omega_{k}}}$ can be implemented using the computation of $\mathbf{S_{ij,\Omega_{k}}}$.

4.3 Application of the proposed method to volume meshes

A natural extension of the proposed method is to apply it to volume meshes. In this case, each tetrahedron can be considered as a unique domain with a given conductivity tensor. In regions where the conductivity is isotropic, tetrahedrons with equal conductivities can be merged. This latter property, following from the fact that the formulation we are proposing is a surfacic one in nature, will warranty substantial savings in terms of the number of unknowns with respect to a standard volumetric method. The surface of the isotropic and homogeneous region is then given by the triangles separating tetrahedron with different conductivity tensors. In practice, if a volume mesh is considered, this requires an additional preprocessing step to extract the interfaces between regions with different conductivity tensors so that the considered geometry is a surfacic mesh. This can be done by comparing for each tetrahedron its conductivity tensor with the conductivity tensor of its neighbor. The two tetrahedrons are merged if the difference between their conductivity tensor is smaller than a chosen threshold $\epsilon$. For example, consider the two tetrahedrons $T_{i}$ and $T_{j}$ with conductivity tensors $\bar{\bar{\sigma}}_{i}$ and $\bar{\bar{\sigma}}_{j}$ respectively, and boundaries given by the set of triangles $\bigcup_{k}t_{ik}$ and $\bigcup_{k}t_{jk}$ respectively. The interface $t_{ij}=\Gamma_{ij}$ between $T_{i}$ and $T_{j}$ is taken into account only if $||\bar{\bar{\sigma}}_{i}-\bar{\bar{\sigma}}_{j}||>\epsilon$, otherwise $T_{i}\cup T_{j}$ will form a unique domain $\Omega_{i}$ such that $\partial\Omega_{i}=\{\bigcup_{k}t_{ik}\cup\bigcup_{k}t_{jk}\}\setminus t_{ij}$. To treat the full volume mesh, the process is repeated for all the tetrahedrons. In the case where the mesh is fully inhomogeneous or when no tetrahedron at all are merged, the system matrix $\mathbf{Z}$ will be extremely sparse since for one line, only seven columns of the matrix are not zero. The reader should keep in mind that using volume mesh is not mandatory to use the presented formulation. It is only a convenient pre-processing step in our situation since, moreover, comparison with volume methods (FEM) will be carried out. Indeed, this step is independent from the integral formulation itself and a pure surface mesh can directly be employed. However, the correct modelling of anisotropic compartments requires to subdivide them. As a consequence, while isotropic compartments are modelled with surfaces, anisotropic ones are virtually using volume elements.

5 Numerical Results

Two sets of numerical experiments have been carried out. Since in case of isotropic conductivity, the presented method is the same as the standard symmetric formulation, we focus on anisotropic conductivity profiles. The first set aimed at assessing the accuracy of the proposed method with canonical anisotropic head meshes (layered anisotropic spheres) while the second set aimed at showing its ability to handle realistic anisotropic conductivity profiles and meshes.

5.1 Modelling of the anisotropy with volume meshes

In modelling the skull layer different approaches are found in the literature. The work of [48] has shown that the different skull layers, a cancellous bone in between two layers of compact bone, do not have the same conductivities. As a consequence, different possibilities are available to model this layer: one can model the skull with a simple isotropic layer, a single anisotropic layer or three layers with different isotropic conductivities. In [49], the authors investigate how to model the skull layer. Their recommendation depends on the available skull geometry accuracy. Obtaining a precise mesh to make out the three layers is not without computational overload and it might explains why modeling the skull with a single anisotropic layer is still very common as for example in [27], [28] or [11]. In this paragraph we present results on spherical layered meshes where the skull is modelled with a single anisotropic layer. The method, however, can of course be applied also in the case of simulations where the skull layer is subdivided into its three isotropic composing layers. The first experiments are based on a three layers sphere mesh of radii 0.87, 0.92 and 1. These layers stand for the brain, skull and scalp respectively. In the first and last layer, the conductivity is set to $1$ while the conductivity in the middle layer is assumed to be radially anisotropic. The conductivity ratio between the first and the second layer is $\sigma_{skull,iso}=1/15$ , following [50]. The ratio between the radial and tangential conductivity in the skull is $1/10$. The radial and tangential conductivity tensor component in each tetrahedra is computed using [27]

[TABLE]

We begin with a volume mesh and extract a surface mesh by comparing the conductivity values in each tetrahedron following the process described in paragraph 4.3. In this fashion, the tetrahedra in the brain and scalp layers are gathered to form only two homogeneous and isotropic domains while the tetrahedra of the skull layer form different domains depending on the accepted conductivity error $\epsilon$ between two neighboring tetrahedra. Figure 4 shows the volume mesh and the surface mesh obtained.

5.1.1 Convergence of the solution

In this first sets of experiment, $\epsilon$ was chosen small enough so that no tetrahedron of the skull layer is merged. This first experiment focused on verifying the convergence of the newly proposed formulation when the number of unknowns is increased. The accuracy is assessed by calculating the relative error of the solution given by the proposed integral formulation with respect to the analytic reference solution available in [13]. The mesh refinement parameter, denoted with $h$, was decreased from $0.38$ to $0.08$ so that the initial number of tetrahedra grew from $3607$ to $112827$. Applying the preprocessing described in 4.3 gave a number of domain increasing from $1495$ to $33117$. A single dipole source with unitary moment along the z-axis was placed in the center of the innermost sphere. The dipole position influence on the numerical results is studied in a second scenario. The results for this test are shown in Figure 5. It is clear that the relative error decreases, showing the convergence and the consistency of the proposed approach. The time that our prototype code needed to solve the anisotropic EEG forward problem in the presented canonical scenario is presented Table 1.

5.1.2 Accuracy for randomized dipole positions and different dipole eccentricities

Since one dipole may not be representative of the complex brain electric sources pattern of activation, we carried out another experiment, with the same anisotropic conductivity parameters, in which dipole positions are randomly set. For each eccentricity, we computed the relative error in 200 simulations of 1 dipole source, for radial and tangential unitary dipole moment. The obtained results, shown in Figure 6, confirm that the proposed approach provides a formulation which is able to handle anisotropic conductivity profile for any dipole position.

5.2 Application to a realistic head mesh

The second experiment is the comparison of the potential obtained in a realistic case with the proposed method and the FEM formulation. Indeed, from grey scale data given by MRI, it is possible to extract volume meshes, (see for example [51]), that are usually used in FEM simulations. From this volume mesh, our method can extract surfaces that delimit domains with the same conductivity. In the presented example, we used [52] to generate a three layers volume head mesh. The fully tetrahedral mesh was used in the FEM solution. A conductivity tensor was associated to each tetrahedra following the choice of the previous anisotropic canonical example. Then, the preprocessing subroutine of our scheme provided a surface mesh for which we applied directly the proposed method. This mesh is shown figure 7(b). The number of tetrahedra was reduced from $187317$ to $51709$ triangular cells. The choice of a merging threshold $\epsilon=0.001$ led to the generation of $11698$ domains. The obtained scalp potential is shown in figure 8. A relative error of $4.5\%$ between the potential computed with the proposed formulation and the fully anisotropic FEM solution was obtained. The potential computed at 256 electrodes positions is presented in figure 9. This figure shows the overlapping curves obtained with the two methods and confirms the validity of the approach. This realistic scenario shows the ability of the proposed method to handle realistic meshes and in particular, if wanted, the skull anisotropy. In case white matter fiber anisotropy is included, no substantial changes must be done: no assumption was made on the geometry of the media. As a consequence, the only difference would reside in the obtention of the mesh. In practice, the number of degrees of freedom may require the use of a fast solver.

6 Discussion and Conclusion

In this paper, a new integral formulation for the EEG forward problem that takes into account anisotropic conductivity profiles was presented. The integral formulation is obtained by leveraging on the representation theorem using an anisotropic fundamental solution. Implementations details are provided and explain the discretization strategy. In particular, a strategy is proposed to tackle inhomogeneous and anisotropic conductivity profile with volume meshes. In this case, the level of inhomogeneity can be controlled by choosing an accepted error in the anisotropic conductivity when building a surfacic mesh from a volume mesh. The level of inhomogeneity does not decrease the accuracy of the solution, on the contrary, it helps to provide a better anatomical model. The numerical results show the correctness of the approach. In the canonical cases, the results given by the new integral equation are matching the analytical solution. In the realistic case, the solution obtained with the proposed formulation is compatible with the one provided by the Finite Elements formulation. One shortcoming of boundary integral formulations, when compared to FEM, is that the arising system matrix is fully populated while in the case of FEM the system matrix is sparse. However, in the presented formulation, as in the symmetric formulation, the number of entries of the system matrix depends on the number of interfaces: the system matrix is block diagonal with the size of the blocks depending on the number of elements in the considered interface. Furthermore, a popular way to fix this disadvantage is to resort to fast solving techniques such as Adaptive Cross Approximation (ACA) or Fast Multipole Method (FMM) that allows to build the system matrix in a linear complexity.

Acknowledgement

This work was supported in part by the CominLabs excellence laboratory, managed by the National Research Agency in the “Investing for the Future” program under reference ANR-10-LABX-07-01, project SABRE and in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC project 321, grant No. 724846).