On the modeling of brain fibers in the EEG forward problem via a new family of wire integral equations

TL;DR

This paper introduces a novel boundary element method (BEM) for EEG forward modeling that accurately incorporates white matter anisotropy, improving the precision of source localization in neuroimaging.

Contribution

A new BEM scheme that models white matter anisotropy efficiently using one-dimensional basis functions, extending traditional isotropic models.

Findings

The new scheme accurately models white matter anisotropy.

It outperforms existing formulations in canonical and realistic cases.

The method is computationally efficient due to boundary formulation.

Abstract

Source localization based on electroencephalography (EEG) has become a widely used neuroimagining technique. However its precision has been shown to be very dependent on how accurately the brain, head and scalp can be electrically modeled within the so-called forward problem. The construction of this model is traditionally performed by leveraging Finite Element or Boundary Element Methods (FEM or BEM). Even though the latter is more computationally efficient thanks to the smaller interaction matrices it yields and near-linear solvers, it has traditionally been used on simpler models than the former. Indeed, while FEM models taking into account the different media anisotropies are widely available, BEM models have been limited to isotropic, piecewise homogeneous models. In this work we introduce a new BEM scheme taking into account the anisotropies of the white matter. The boundary nature of the formulation allows for an efficient discretization and modelling of the fibrous nature of the white matter as one-dimensional basis functions, limiting the computational impact of their modelling. We compare our scheme against widely used formulations and establish its correctness in both canonical and realistic cases.

Figures

Tables

Table 1. Table 1 : Definition of the different coefficients associated with the discretization of the symmetric formulation.

Condition	${\alpha }_{ij}$	${\beta }_{ij}$	${\gamma }_{ij}$	${\theta }_j$
$j=1$	$-$	$-$	$-$	1
$j=i$	${\sigma }_i+{\sigma }_j$	$-2$	${\sigma }_i^{-1}+{\sigma }_j^{-1}$	$-$
$j=i-1$	$-{\sigma }_i$	$1$	$-{\sigma }_i^{-1}$	$-$
$j=i+1$	$-{\sigma }_j$	$1$	$-{\sigma }_j^{-1}$	$-$
otherwise	$0$	$0$	0	0

Table 2. Table 2 : Coordinates of the wires in the x ​ y 𝑥 𝑦 xy plane corresponding to Figure   6 .

x	-0.4	-0.2	0	0.2	0.4	-0.2	0	0.2	-0.2	0	0.2
y	0	0	0	0	0	-0.2	-0.2	-0.2	0.2	0.2	0.2

Equations

$$\nabla \cdot \left(\underline{\underline{\sigma }}(\mathbf{r})\nabla \varphi (\mathbf{r})\right)=\nabla \cdot \mathbf{J}(\mathbf{r}),\mathbf{r}\in \Omega ,$$

$$\nabla \cdot \left(\underline{\underline{\sigma }}(\mathbf{r})\nabla \varphi (\mathbf{r})\right)$$

$${\sigma }_i\text{upDelta}\varphi (\mathbf{r})$$

$${\left[\varphi \right]}_j$$

$${\left[\sigma {\partial }_{\hat{\mathbf{n}}}\varphi \right]}_j$$

$$\mathbf{J}(\mathbf{r})=\mathbf{P}\delta (\mathbf{r}-{\mathbf{r}}_0),$$

$$v_{\text{dip}}(\mathbf{r})=\frac{\mathbf{P}\cdot (\mathbf{r}-{\mathbf{r}}_0)}{4\pi {\sigma }_1{\left|\mathbf{r}-{\mathbf{r}}_0\right|}^3}.$$

$$\text{upDelta}\varphi (\mathbf{r})=\nabla \cdot \left(\frac{\mathbf{J}(\mathbf{r})}{{\sigma }_1}+\underline{\underline{\kappa }}(\mathbf{r}){\mathbb{J}}_f(\mathbf{r})\right),\mathbf{r}\in {\Omega }_1,$$

$$\underline{\underline{\kappa }}\left(\mathbf{r}\right)=\left({\underline{\underline{\sigma }}}^{-1}(\mathbf{r})-{{\sigma }_1}^{-1}\underline{\underline{I}}\right),$$

$$\underline{\underline{\kappa }}(\mathbf{r})=({\sigma }_f^{-1}-{\sigma }_1^{-1})\hat{\mathbf{t}}(\mathbf{r}){\hat{\mathbf{t}}}^{\mathrm{T}}(\mathbf{r})$$

$$(\mathcal{S}{u|}_{{\Gamma }_i})(\mathbf{r})={\int }_{{\Gamma }_i}G(\mathbf{r}-{\mathbf{r}}^{\prime })u({\mathbf{r}}^{\prime })\mathrm{d}s({\mathbf{r}}^{\prime }),\mathbf{r}\in \Omega ,$$

$$(\mathcal{D}{u|}_{{\Gamma }_i})(\mathbf{r})={\int }_{{\Gamma }_i}{\partial }_{{\hat{\mathbf{n}}}^{\prime }}G(\mathbf{r}-{\mathbf{r}}^{\prime })u({\mathbf{r}}^{\prime })\mathrm{d}s({\mathbf{r}}^{\prime }),\mathbf{r}\in \Omega ,$$

$$(\mathcal{K}{u|}_{{\Gamma }_i})(\mathbf{r})={\int }_{{\Gamma }_i}{\partial }_{\hat{\mathbf{n}}}G(\mathbf{r}-{\mathbf{r}}^{\prime })u({\mathbf{r}}^{\prime })\mathrm{d}s({\mathbf{r}}^{\prime }),\mathbf{r}\in {\Gamma }_j,$$

$$(\mathcal{N}{u|}_{{\Gamma }_i})(\mathbf{r})={\int }_{{\Gamma }_i}{\partial }_{\hat{\mathbf{n}}}{\partial }_{{\hat{\mathbf{n}}}^{\prime }}G(\mathbf{r}-{\mathbf{r}}^{\prime })u(\mathbf{r})\mathrm{d}s({\mathbf{r}}^{\prime }),\mathbf{r}\in {\Gamma }_j,$$

$$G(\mathbf{r}-{\mathbf{r}}^{\prime })=\frac{1}{4\pi \left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}$$

$$(\mathcal{V}\mathbf{u})(\mathbf{r})={\int }_{\Lambda }g(\mathbf{r}-{\mathbf{r}}^{\prime })\nabla \cdot (\underline{\underline{\kappa }}\mathbf{u}({\mathbf{r}}^{\prime }))\mathrm{d}l({\mathbf{r}}^{\prime }),\mathbf{r}\in \Omega ,$$

$$g(\mathbf{r}-{\mathbf{r}}^{\prime })={\int }_0^a{\int }_0^{2\pi }G(\mathbf{r}-{\mathbf{r}}^{\prime }){\rho }^{\prime }\mathrm{d}{\rho }^{\prime }\mathrm{d}{\theta }^{\prime },$$

$$(\mathcal{D}{\varphi |}_{\partial {\Omega }_i})(\mathbf{r})-(\mathcal{S}{\xi |}_{\partial {\Omega }_i})(\mathbf{r})={\int }_{{\Omega }_i}\varphi (\mathbf{r})\text{upDelta}G(\mathbf{r}-{\mathbf{r}}^{\prime })-G(\mathbf{r}-{\mathbf{r}}^{\prime })\text{upDelta}\varphi (\mathbf{r})\mathbf{r}\in {\Omega }_i\setminus {\Gamma }_i$$

$$(\mathcal{D}{\left.\kern-1.2pt\phi\vphantom{|}\right|_{\partial{\Omega}_{i}}})(\bm{r})-(\mathcal{S}{\left.\kern-1.2pt\xi\vphantom{|}\right|_{\partial{\Omega}_{i}}^{\pm}})(\bm{r})=\\ \phi(\bm{r})+\left\{\begin{array}[]{cl}v_{\text{dip}}(\bm{r})-\sum\limits_{k=1}^{{N_{\text{f}}}}{\mathcal{V}\bm{\mathbb{J}}_{\Lambda_{\text{k}}}}(\bm{r})&\bm{r}\in{\Omega_{i}}\setminus\Gamma_{i},\,i=1\,,\\ 0&\bm{r}\in{\Omega_{i}}\setminus\Gamma_{i},\,i=2\dots N\,,\end{array}\right.$$

$$(\mathcal{D}{\left.\kern-1.2pt\phi\vphantom{|}\right|_{\partial{\Omega}_{i}}})(\bm{r})-(\mathcal{S}{\left.\kern-1.2pt\xi\vphantom{|}\right|_{\partial{\Omega}_{i}}^{\pm}})(\bm{r})=\\ -\frac{1}{2}{\left.\kern-1.2pt\phi\vphantom{|}\right|_{\partial{\Omega}_{i}}}(\bm{r})+\left\{\begin{array}[]{cl}v_{\text{dip}}(\bm{r})-\sum\limits_{k=1}^{{N_{\text{f}}}}{\mathcal{V}\bm{\mathbb{J}}_{\Lambda_{\text{k}}}}(\bm{r})&\bm{r}\in{\partial\Omega_{i}},\,i=1\,,\\ 0&\bm{r}\in{\partial\Omega_{i}},\,i=2\dots N\,,\end{array}\right.$$

$$(\mathcal{N}{\left.\kern-1.2pt\phi\vphantom{|}\right|_{\partial{\Omega}_{i}}})(\bm{r})-(\mathcal{K}{\left.\kern-1.2pt\xi\vphantom{|}\right|_{\partial{\Omega}_{i}}^{\pm}})(\bm{r})=\\ -\frac{1}{2}{\left.\kern-1.2pt\xi\vphantom{|}\right|_{\partial{\Omega}_{i}}^{\pm}}(\bm{r})+\left\{\begin{array}[]{cl}v_{\text{s}}(\bm{r})-\sum\limits_{k=1}^{{N_{\text{f}}}}{\mathcal{W}\bm{\mathbb{J}}_{\Lambda_{\text{k}}}}(\bm{r})&\bm{r}\in{\partial\Omega_{i}},\,i=1\,,\\ 0&\bm{r}\in{\partial\Omega_{i}},\,i=2\dots N\,,\end{array}\right.$$

$$(\mathcal{W}\bm{u})(\bm{r})=\int_{\Lambda}\partial_{\hat{\bm{n}}}g(\bm{r}-\bm{r}^{\prime})\nabla\cdot\left(\underline{\underline{\kappa}}\bm{u}(\bm{r}^{\prime})\right)\,\mathrm{d}l(\bm{r}^{\prime})\$$

$$v_{\text{s}}(\mathbf{r})={\partial }_{\hat{\mathbf{n}}}{v}_{\text{dip}}(\mathbf{r}).$$

$$\begin{gathered} {\sigma }_1{v}_{\text{dip}}(\mathbf{r})-{\sigma }_1\sum _{k=1}^{N_{\text{f}}}\mathcal{V}{\mathbb{J}}_{{\Lambda }_{\text{k}}}(\mathbf{r})=\frac{{\sigma }_j+{\sigma }_{j+1}}{2}\varphi (\mathbf{r})-\sum _{i=1}^N\left({\sigma }_{i+1}-{\sigma }_i\right)(\mathcal{D}{\varphi |}_{{\Gamma }_i})(\mathbf{r}), \\ \mathbf{r}\in {\Gamma }_j,j=1\dots N, \end{gathered}$$

$$\begin{gathered} {\sigma }_1\nabla v_{\text{dip}}(\mathbf{r})-{\sigma }_1\nabla \sum _{k=1}^{N_{\text{f}}}\mathcal{V}{\mathbb{J}}_{{\Lambda }_{\text{k}}}(\mathbf{r})={\sigma }_1{\underline{\underline{\sigma }}}_{{\Lambda }_n}^{-1}{\mathbb{J}}_{{\Lambda }_{\text{n}}}(\mathbf{r})-\nabla \sum _{i=1}^N\left({\sigma }_{i+1}-{\sigma }_i\right)(\mathcal{D}{\varphi |}_{{\Gamma }_i})(\mathbf{r}), \\ \mathbf{r}\in {\Lambda }_n,n=1\dots N_{\text{f}}. \end{gathered}$$

$$v_{\text{s}}(\mathbf{r})-\sum _{k=1}^{N_{\text{f}}}\mathcal{W}{\mathbb{J}}_{{\Lambda }_{\text{k}}}(\mathbf{r})=\frac{1}{2}({\xi |}_{{\Gamma }_j}^{-}+{\xi |}_{{\Gamma }_j}^{+})(\mathbf{r})-\sum _{i=1}^N\mathcal{K}({\xi |}_{{\Gamma }_i}^{-}-{\xi |}_{{\Gamma }_i}^{+})(\mathbf{r}),\mathbf{r}\in {\Gamma }_j,j=1\dots N,$$

$$q_{{\Gamma }_j}={\xi |}_{{\Gamma }_j}^{-}-{\xi |}_{{\Gamma }_j}^{+}=\left(\frac{{\sigma }_{j+1}-{\sigma }_j}{{\sigma }_{j+1}}\right){\xi |}_{{\Gamma }_j}^{-},$$

$${\xi |}_{{\Gamma }_j}^{-}+{\xi |}_{{\Gamma }_j}^{+}=\left(\frac{{\sigma }_{j+1}+{\sigma }_j}{{\sigma }_{j+1}-{\sigma }_j}\right){q|}_{{\Gamma }_j}.$$

$$v_{\text{s}}(\mathbf{r})-\sum _{k=1}^{N_{\text{f}}}\mathcal{W}{\mathbb{J}}_{{\Lambda }_{\text{k}}}(\mathbf{r})=\frac{{\sigma }_j+{\sigma }_{j+1}}{2({\sigma }_{j+1}-{\sigma }_j)}{q|}_{{\Gamma }_j}(\mathbf{r})-\sum _{i=1}^N(\mathcal{K}{q|}_{{\Gamma }_i})(\mathbf{r}),\mathbf{r}\in {\Gamma }_j,j=1\dots N.$$

$$\varphi (\mathbf{r})=v_{\text{dip}}(\mathbf{r})-\sum _{k=1}^{N_{\text{f}}}\mathcal{V}{\mathbb{J}}_{{\Lambda }_{\text{k}}}(\mathbf{r})+\sum _{i=1}^N(\mathcal{S}{q|}_{{\Gamma }_i})(\mathbf{r}),$$

$$\nabla v_{\text{dip}}(\mathbf{r})-\nabla \sum _{k=1}^{N_{\text{f}}}\mathcal{V}{\mathbb{J}}_{{\Lambda }_{\text{k}}}(\mathbf{r})={\underline{\underline{\sigma }}}_{{\Lambda }_n}^{-1}{\mathbb{J}}_{{\Lambda }_{\text{n}}}(\mathbf{r})-\nabla \sum _{i=1}^N(\mathcal{S}{q|}_{{\Gamma }_i})(\mathbf{r}),\mathbf{r}\in {\Lambda }_{\text{n}},n=1\dots N_{\text{f}}.$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Electroencephalography (EEG) based source localization has gained an increasing popularity as a reliable neuroimaging modality in research and medical practice [1, 2, 3]. Using scalp measured potentials, various algorithms have been proposed for the retrieval of the location of the neuro-generators [4]. Many of these algorithms rely on an accurate solution of the associated forward problem which maps a given setting of sources and head model to the corresponding scalp potential. The complexity of the head geometry and its underlying conductivity, however, precludes the use of analytical methods and one has to adopt numerical approximations. With their renowned high accuracy and robustness, integral equations-based methods remain the preferred choice for researchers [5, 6]. In particular, the boundary element method (BEM) only requires the discretization of the boundaries, thus reducing the overall dimensionality. Moreover, given the smoothness of its underlying kernel, it is possible to augment BEM with fast algorithms such as the adaptive cross approximation (ACA) or the fast multipole method (FMM) [7, 8], which further reduce its computational complexity. The three most widely employed BEM formulations for the EEG forward problem are the adjoint double layer (ADL), the double layer (DL) and the symmetric (SY) approaches [9, 10, 11]. By leveraging on methods of layer potentials, these methods solve Poisson equation under the assumption of isotropic media [12]. The DL formulation is a direct approach in which the potential is obtained directly while the ADL formulation solves first for an auxiliary unknown before integrating it to obtain the electric potential. Differently from the two previous approaches, the SY formulation simultaneously involves two surface unknowns. Despite its larger system of equations, it has a block diagonal structure [12]. For more details of these methods, their relative merits and disadvantages, the reader is referred to [13, 10, 11].

Despite their advantages, BEM-based formulations are restricted to isotropic and piece-wise homogeneous problems. This is a significant limitation since white matter anisotropy has a considerable impact [14, 15] on the accuracy of source localization procedures. These early results have been obtained with differential based methods and entire volume discretization [14, 15], which is computationally expensive. More recently, integral techniques accounting for the white matter anisotropy have been introduced [16, 13]; they do however also require discretization of the entire head volume.

The anisotropy of the white matter tissue arises from its underlying assembly of bundles of parallelly-oriented axon [17, 18]. This suggests that the apparent inhomogeneous anisotropy is actually structured and may be expressed in terms of these axons’ fibers. This observation has been leveraged on in [19] by replacing a single fiber by dipolar sources of constant magnitude. The forward problem was subsequently solved iteratively with the symmetric formulation. However, this work does not account for the coupling and interactions between different fibers which is essential for precise forward solution.

The work presented in this paper aims at extending the three main BEM (EEG) formulations to take into account the anisotropic and inhomogeneous conductivity of the white matter. This is achieved by a modelization of the white matter connectivity. Indeed, using diffusion weighted MRI (DW-MRI) it is possible to track axon fibers and reveal the underlying network of the white matter [20]. One-dimensional basis functions are used for the modelization of the fibers which results in efficient and accurate forward solutions. As a byproduct, the new technique we present could further improve the recently introduced approaches exploiting the brain connectivity patterns in source estimation [21, 22]. Some preliminary results have been presented in [23]. Several numerical experiments validate the new schemes in canonical and realistic settings.

The reader should note that 1D formulations have been extensively studied in the context of high frequency electromagnetic modeling of wire-like structures [24, 25, 26], although those schemes, for perfect electrically conducting wires, are only mildly related to the ones presented here.

The paper is organized as follow: the notations is set and some background is recalled in Section 2; the new equations and their discretizations are then derived in Section 3 and Section 4, respectively. The new schemes are validated with various simulations and tests in Section 5 before closing with conclusions in Section 6.

2 Background and notations

Consider an electric volume current density $\bm{J}$ residing in a conducting medium $\Omega\subset\mathbb{R}^{3}$ composed of $N$ nested, piecewise homogeneous sub-regions $\Omega_{\text{i}}$ such that $\Omega=\bigcup\nolimits_{i=1}^{N}{{\Omega_{\text{i}}}}$ with $\Omega_{1}$ being the innermost layer. Each sub-region is associated with an isotropic conductivity $\sigma_{i}$ and delimited by the Lipschitz surface $\Gamma_{i}$ with $\Gamma=\bigcup\nolimits_{i=1}^{N}{{\Gamma_{\text{i}}}}$ and $\Gamma_{i}\cup\Gamma_{i-1}=\partial\Omega_{i}$ (Figure 1).

Let $\Lambda_{\text{j}}$ be a curve modeling a bundle of parallel white matter fibers and $\Lambda=\bigcup_{j=1}^{N_{\text{f}}}{{\Lambda_{\text{j}}}}$ their union. These fibers assume a thin cylindrical shape of circular section $a$ [27, 28] and potentially contain junctions; they give rise to a tensorial conductivity profile in which the conductivity along the fiber $\sigma_{f}$ is different from the conductivity in the transversal direction $\sigma_{1}$. In the quasi-static regime, the electric potential $\phi$ is related to the current density $\bm{J}$ via Poisson’s equation

[TABLE]

where the local electric conductivity $\underline{\underline{\sigma}}$ is described by a $3\times 3$ symmetric tensor [29]. For the sake of simplicity in explaining our new method, we neglect the anisotropy of the skull and thus, given that the neuron’s fibers are only present in the innermost region $\Omega_{1}$, eq. 1 can be rewritten as

[TABLE]

where it was assumed that the current sources are present only in the innermost layer that corresponds to the brain. Two transmission conditions are associated with each of these equations: (i) a Dirichlet condition that enforces the continuity of the electric potential across interfaces and (ii) a Neumann condition that enforces the continuity of the electric current flux, i.e.

[TABLE]

where the bracket notation $[g]_{j}$ denotes the jump of a function $g$ across $\Gamma_{j}$ and $\partial_{\hat{\bm{n}}}g=\hat{\bm{n}}\cdot\nabla g$ with $\hat{\bm{n}}=\hat{\bm{n}}(\bm{r})$ is the unit vector normal to $\Gamma_{j}$ pointing outward of $\Omega_{j}$. Note that the fibers do not come into contact with the inner surface. In the context of the EEG forward problem, brain sources are commonly modeled as current dipoles [30, 31], i.e.

[TABLE]

in which $\bm{P}$ and $\delta$ respectively denote the dipole moment and the Dirac delta function. The electric potential induced by this elementary source in an infinite homogeneous domain of conductivity $\sigma_{1}$ reads

[TABLE]

3 Integral equation based formulations

By transforming Poisson equation into an integral equation, conventional BEM formulations (SL, DL and SY), have been particularly attractive as they offer computational savings in comparison with other alternatives. In general, reformulating a partial differential equation as an integral expression requires knowledge of its fundamental solution. When considering the anisotropy of the white matter however, eq. 2 involves position and orientation dependent tensors for which the corresponding fundamental solution does not exist in closed form, for general geometries. Equation 2 should then be recast into an equivalent one at a reduced dimensionality by extracting the Laplacian operators and using Green’s identities. In particular, this choice not only allows for a unified treatment of eqs. 2 and 3, but also reduces the effect of the anisotropy to an one-dimensional apparent volume current density along the fibers. Consequently, the framework of standard BEM formulations can be extended to handle the anisotropy of the white matter. To that end, eq. 2 can be re-expressed as

[TABLE]

where $\bm{\mathbb{J}}_{f}=\underline{\underline{\sigma}}\nabla\phi$ is the apparent volume current density along the fibers and $\underline{\underline{\kappa}}$ is the conductivity contrast defined as

[TABLE]

in which $\underline{\underline{I}}$ is the identity tensor. Note that $\underline{\underline{\kappa}}$ is zero everywhere except on the fibers where it has the form

[TABLE]

with $\hat{\bm{t}}(\bm{r})$ being the unit vector tangential to $\Lambda_{j}$. We remind the reader that it was assumed that the conductivity of the fibers is ${\sigma_{1}}$ transversally and $\sigma_{f}$ longitudinally. In order to derive an integral representation for the potential and its flux using Green’s second identity, the following well known operators are introduced:

the single layer operator

[TABLE]

• 2.

the double layer operator

[TABLE]

• 3.

the adjoint double layer operator

[TABLE]

• 4.

the hypersingular operator

[TABLE]

where

[TABLE]

is the fundamental solution associated with the Laplacian. In addition, a new operator is introduced to handle the fiber contributions

[TABLE]

where the associated wire kernel is defined as

[TABLE]

and $\rho$ and $\theta$ are the usual polar coordinates in the fiber’s transverse plane.

The starting point of our development is Green’s second identity, which states that:

[TABLE]

where $\xi=\partial_{\hat{\bm{n}}}\phi$ is the derivative of the potential in the normal direction. Using eqs. 3 and 8 and the property of the fundamental solution, eq. 18 reduces to the following

[TABLE]

Taking the limit $\bm{r}\rightarrow\partial\Omega$, the following integral representation for the electric potential is derived

[TABLE]

By differentiating eq. 19 with respect to $\bm{r}$ in the direction normal to the boundary, an integral representation for the potential flux can be obtained

[TABLE]

where

[TABLE]

and

[TABLE]

It is worth noting that eqs. 20 and 21 are written for a normal vector pointing outward. A consistent change of signs should be made when the normals are pointing inward, which is the case for $\Gamma_{i-1}$. In the inner most layer $\Omega_{1}$ the last term of the right-hand side in eqs. 20 and 21 represents the effect of the fibers; it describes the local anisotropic conductivity.

3.1 Double layer-Wire formulation

Equations 20 and 21 have two surface unknowns, one of which could be discarded by using the boundary conditions eqs. 4 and 5; depending on the variable discarded two different formulations can be obtained. The double layer-wire formulation is obtained if the surface electric potential $\phi(\bm{r})$ is the remaining unknown. This formulation can be derived after multiplying eq. 20 with the local conductivity and summing the contribution of all the regions $\Omega_{i}$

[TABLE]

where the $\mathcal{S}$ operator term cancels out by enforcing the transmission condition (5).

Equation 24 simultaneously involves the surface potential $\phi$ and the current density $\bm{\mathbb{J}}$ as unknowns, and therefore needs to be complemented with a second equation. The second equation is obtained by applying the gradient operator to eq. 24

[TABLE]

Combining eqs. 25 and 24 constitutes the first new formulation and will be referred to as the double layer-wire (DLW) formulation.

3.2 Single layer-Wire formulation

Differently from the DLW that is formulated in terms of surface potentials, the single layer-Wire (SLW) formulation is derived from eq. 21 and solves for the jump of the potential’s normal derivative across an interface. Thus, applying eq. 21 to each region $\Omega_{j}$ and summing up their contributions yields

[TABLE]

where the $\mathcal{N}$ operator term cancels out by enforcing the transmission condition (4). After introducing

[TABLE]

the difference between normal derivatives can be expressed as

[TABLE]

Substituting back eq. 28 in eq. 26 forms the single layer formulation

[TABLE]

Similarly to the DLW eq. 29 exhibits two unknowns and needs to be complemented. The complementary equation will be derived from eq. 20 by summing the contributions of all regions

[TABLE]

where $\mathcal{D}$ vanishes due to condition eq. 4. A current equation is obtained by applying the gradient operator to eq. 30

[TABLE]

Subsequently to finding $q$, the electric potential can be computed via eq. 30.

3.3 Symmetric-Wire formulation

The symmetric formulation leverages on a combination of eq. 20 and eq. 21 applied, in contrast with the two previous formulations, to adjacent regions only. Summing these contributions yields

[TABLE]

The current flux ${\left.\kern-1.2ptd\vphantom{|}\right|_{\Gamma_{i}}}=\sigma_{i}\partial_{\hat{\bm{n}}}{\left.\kern-1.2pt\xi\vphantom{|}\right|_{\Gamma_{i}}^{-}}=\sigma_{i+1}\partial_{\hat{\bm{n}}}{\left.\kern-1.2pt\xi\vphantom{|}\right|_{\Gamma_{i}}^{+}}$ (by virtue of condition eq. 5), can be substituted in eq. 32

[TABLE]

This expression constitutes the first equation of the symmetric formulation. It has three unknowns, the surface potential, the normal component of the surface current density and the fibers current density. Therefore, two other equations are needed. In order to derive a second equation, eq. 33 is multiplied by the local conductivity and applied to adjacent domains, yielding

[TABLE]

For the third equation, the gradient operator of eq. 20 is applied to the innermost layer, which leads to the current equation

[TABLE]

Note that the quantities restricted to non-existing surfaces i.e. $\Gamma_{0}$ and $\Gamma_{N+1}$ are set to zero and that on the outermost layer ${\left.\kern-1.2ptd\vphantom{|}\right|_{\Gamma_{N}}}$ is identically Zero. This formulation requires the solution of two surface equations out of which the surface unknowns interact with only their immediate neighbors. This will give rise to a block diagonal matrix, thus reducing the apparently higher computational cost.

4 Discretization

The numerical solution of the presented equations is achieved following a Galerkin approach. In this respect, the different head surfaces $\Gamma_{i}$ are tessellated into triangular meshes and the fibers $\Lambda_{j}$ into cylindrical segments. On these finite elements, each unknown $S(\bm{r})$ is approximated by a linear combination of the $N_{x}$ basis functions $\{x_{i}\}$

[TABLE]

where $a_{i}=\left<S(\bm{r}),x_{i}(\bm{r})\right>$. In order to obtain a square linear system, the discretized equations are then tested with an appropriate set of functions of same cardinality as the set of basis functions. The choice of these finite elements is not arbitrary and must be in accordance with the operators’ mapping properties i.e. the basis functions should span the domain of the operator and the testing functions should span the dual of its range [11, 32]. The functions used to discretize the different unknowns must be capable of satisfying their different physical properties, for instance the discretization of the current density should not permit the existence of jumps. In this paper we considered patch $\{\varphi_{n}(\bm{r})\}$ and pyramid $\{\psi_{n}(\bm{r})\}$ functions to expand the surface unknowns $\phi$, $q$ and $d$ depending of the formulation and hat basis functions $\{\bm{\lambda}_{n}(\bm{r})\}$ to expand the current density $\mathbb{J}(\bm{r})$. The pyramid and patch basis functions $\psi_{n}(\bm{r})$ are respectively expressed as

[TABLE]

where $\bm{r_{n}},\bm{r_{i}},\bm{r_{j}}$ are the position vectors of the vertices constituting the triangle $T\text{r}_{n}$. Figures 2(b) and 2(a) presents the schematic definitions of these basis functions.

The current density is expanded with oriented hat functions whose support are the cylindrical segments $s_{k}=(\bm{r_{k}};\,\bm{r_{k+1}})$ and $s_{k+1}=(\bm{r_{k+1}};\,\bm{r_{k+2}})$ (Figure 3) and defined as

[TABLE]

It should be noted that the hat basis functions are continuous and thus automatically enforce the jump condition of the current density.

4.1 Discretization of the double layer-wire formulation

In eqs. 24 and 25, the surface potential $\phi$ is discretized with pyramid basis functions and the current density $\bm{\mathbb{J}}$ is discretized with hat basis functions. Equations 24 and 25 are then tested with pyramid and hat functions respectively. This gives rise to the following matrix system

[TABLE]

where the matrix entries are

[TABLE]

and, ${\left\langle{f,g}\right\rangle_{x}}=\int_{x}{f\cdot g}\,\,\mathrm{d}x$ denotes the duality product. The entries of the right-hand side are

[TABLE]

4.2 Discretization of the single layer-wire formulation

Similarly to the previous approach, the surface unknown $q$ in eqs. 29 and 30 is discretized with pyramid basis functions and the current density $\bm{\mathbb{J}}$ is discretized with hat basis functions. Equations 29 and 30 are then tested with pyramid and hat functions respectively. This gives rise to the system

[TABLE]

where the matrix entries are

[TABLE]

and where the entries of the right-hand side are

[TABLE]

4.3 Discretization of the symmetric-wire formulation

In contrast with the two previous approaches, the symmetric formulation (eqs. 33, 34 and 35) has two surface unknowns: the potential $\phi$ which is discretized with pyramid basis functions, the current flux with patch basis functions and the current density $\bm{\mathbb{J}}$ with hat basis functions. Equations 33, 34 and 35 are tested with patch, pyramid and hat basis functions respectively, resulting in the following system of equations

[TABLE]

in which the system entries are defined as

[TABLE]

the entries of the right-hand side are

[TABLE]

and the coefficients $\alpha$, $\beta$, $\gamma$ and $\theta$ are defined in Table 1.

5 Numerical results

In this section the newly developed integral formulations are validated and their performances are studied through several numerical examples. The parameters of the simulations are given in normalized units.

5.1 Convergence of the solution

In order to demonstrate that the proposed formulations are capable of capturing the anisotropic conductivity caused by the brain fibers and do converge to the exact solution, we have simulated a cubic block whose anisotropic conductivity is $10$ along the $z$ axis and $1$ in the $(x,\,y)$ plane, residing inside a three layered sphere (Figure 4(a)). The radii of the spheres are 0.87$,$0.92 and $1$ respectively. The cube, whose side length is equal to $0.7$, is placed at their center. The conductivity of the different spheres are $1$, $1\text{/}15$ and $1$ respectively. A current source with a dipolar moment equal to $[$1$,$1$,$1$]$ is set at $[$0.4$,$0$,$0$]$. In order to account for the anisotropic effect of the cube with our formulations, we have created a grid of $64$ equally spaced fiber, as illustrated in fig. 4(b). The wires have a radius of $0.05$ and their conductivity is set to be $10$ along the wires and $1$ in their transverse direction. A convergence analysis has been carried out in which the model is discretized with increasingly refined mesh (the number of wires has been kept constant). Note that the exact wire structure shown in Figure 4(b) could have been solved with FEM and used as reference rather than the one obtained with the cube structure in Figure 4(a), however a more extreme case was preferred to illustrate the merits of the new schemes by choosing a different modelization of the underlying physics.

Figure 5 reports the obtained relative error as a function of the mesh edge length; a FEM solution corresponding to highly refined mesh is used as reference. It is clear that the three formulations converge to the reference solution and can indeed account for the anisotropy of the medium.

5.2 Accuracy for different dipole eccentricities

In the second test, we have studied the effect of source eccentricity on the computed potential. Three concentric spheres of radius $0.87$, $0.92$ and $1$ have been considered. Eleven vertical fibers of radius $0.05$ were set at the coordinates summarized in Table 2.

The conductivities of the different layers of the sphere were set to $1$, $1\text{/}15$ and $1$ and the fibers to $10$ along the $z$ direction and $1$ in the transversal direction. The model was discretized with $642$ nodes per surface and $15$ segments per fiber. The forward problem was then solved for a varying dipole position: along and away from the fibers as shown in fig. 6 with red dots. The computed relative error, where a high resolution FEM was used as a reference, is shown in figs. 7(b) and 7(a) for the two cases. In order to illustrate the error introduced when neglecting the anisotropic conductivity of the fibers, we have also included the relative error produced by the analytic solution of the same spherical model in the absence of the fibers.

In general the accuracy of the three numerical solutions decreases for shallow sources. This behavior is due to the singularity of the source and is in agreement with what has been reported in the literature [33, 12, 34]. It is also observed that not accounting for the anisotropic conductivity of the fibers leads to higher errors, especially in vicinity of the fibers. As expected, these errors decrease when the source is moved away from the fibers (fig. 7(b)) and remains stable when moving in their vicinity (fig. 7(a)).

5.3 Application to a realistic head mesh

As a last numerical test, a realistic head model obtained from MRI images is considered (Figure 8(a)). Using standard procedures (see for example [35]), we have constructed a 3 layered mesh in which each domain represents the brain, the skull and the scalp, each of which is made of $6248$, $8328$ and $9346$ triangles, respectively. Furthermore, the white matter fibers are recovered using DTI-based tractography implemented in [36]. The conductivity of the different tissues is set to $0.33$, $0.067$ and $0.33$ for the scalp, skull and brain respectively. The conductivity of the fibers is set to be $0.33$ in their local transverse direction and $10$ times greater in their longitudinal direction. Following the EGI system [37], a set of $256$ electrodes has been placed on the scalp as shown in Figure 8(b). At these positions, the electric potential was computed using the newly introduced schemes. For the sake of comparison, we have also computed the solution with FEM, on a volume mesh of 10 million tetrahedrons. We show the results obtained in Figure 9, where we observe that the four formulations are in agreement. In Figure 8(c), we plot the magnitude of the current density along the fibers.

6 Conclusion

The correct modeling of the electric properties of the head is crucial for an accurate forward solution and, consequently, for brain source reconstruction. This includes the anisotropic behavior of the white matter, given its impact on the scalp potential. In this paper, we have presented new integral techniques that can handle the anisotropic conductivity profile of the head and thus extend the application of conventional BEM approaches. The one dimensional nature of the wire basis functions ensure the computational efficiency of the schemes. It has been shown throughout several numerical tests that the computed potential exhibit high accuracy and stability making it a competitive alternative to differential equations based methods.

Acknowledgment

This work has been funded 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).