# A Hybrid Volume-Surface-Wire Integral Equation for the Anisotropic Forward
Problem in Electroencephalography

Maxime Y. Monin  Student Member, IEEE, Lyes Rahmouni  Adrien Merlini  Member,
IEEE, and Francesco P. Andriulli  Senior Member, IEEE

###### Abstract

Solving the electroencephalography (EEG) forward problem is a fundamental step
in a wide range of applications including biomedical imaging techniques based
on inverse source localization. State-of-the-art electromagnetic solvers
resort to a computationally expensive volumetric discretization of the full
head to account for its complex and heterogeneous electric profile. The more
efficient, popular in biomedical imaging circles, but unfortunately
oversimplifying Boundary Element Method (BEM) relies instead on a piecewise-
uniform approximation that severely curbs its application in high resolution
EEGs. This contribution lifts the standard BEM contraints by treating the
local anisotropies with adequate wire and thin volume integral equations that
are tailored to specific structures of the fibrous white matter and the
inhomogeneous skull. The proposed hybrid integral equation formulation thereby
avoids the full volumetric discretization of the head medium and allows for a
realistic and efficient BEM-like solution of the anisotropic EEG forward
problem. The accuracy and flexibility of the proposed formulation is
demonstrated through numerical experiments involving both canonical and
realistic MRI-based head models.

Anisotropic Forward Problem, Electroencephalography (EEG), Electromagnetic
Integral Equations, Tractography.

22footnotetext: This work was supported by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 724846, project 321). Data were provided in part by the
Human Connectome Project, WU-Minn Consortium (Principal Investigators: David
Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and
Centers that support the NIH Blueprint for Neuroscience Research; and by the
McDonnell Center for Systems Neuroscience at Washington
University.33footnotetext: The authors are with the Department of Electronics
and Telecommunications, Politecnico di Torino, Turin, 10129, Italy (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>francesco.andriulli@polito.it).

©2020 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

##  1 Introduction

High-resolution electroencephalography (EEG) is a staple neuroimaging tool
employed in a wide variety of biomedical applications ranging from clinical
diagnostics to research in neuroscience. An EEG system records the electric
potential on scalp electrodes generated by the brain neural activity at a high
temporal resolution and is both affordable and portable, which has facilitated
its broad adoption in the last few decades [1, 2, 3]. Compared to other
imaging modalities, however, EEG imaging is sensitive to the smearing effects
of various head tissues, resulting in low spatial resolution. This challenge
can be tackled through advanced computational techniques [4], thus enabling
precise analysis of deep neural activity with EEG [5].

The EEG forward problem is a characterization of the relationship between the
neural activity in the brain and the electric potential measured on the scalp.
Its solution is necessary for the associated inverse problem [6] that consists
in inferring knowledge about brain activity from EEG data. This
characterization is challenging due to the complex nature of the human head.
In particular, the skull is made of varying thicknesses of soft bone and
highly resistive hard bone layers [7, 8], whereas diffusion Magnetic Resonance
Imaging (dMRI) shows that the brain white matter has a fiber-like structure of
coherently oriented axon bundles connecting different brain regions [9].
Overall, these fine structures make the head volume strongly inhomogeneous and
anisotropic, and must be suitably taken into account as they have a strong
influence on brain imaging accuracy [10, 11, 12]. In realistic applications,
the forward problem is solved using numerical techniques [13] derived from
full-wave or quasi-static differential and integral formulations of Maxwell’s
equations [14, 15, 16, 17, 18, 19]. State-of-the-art differential methods such
as the Finite Element Method (FEM) can easily incorporate the local variations
of tissue conductivity, but rely on a computationally expensive volumetric
discretization of the entire head. The Boundary Element Method (BEM) [20, 21]
is a popular alternative which reformulates the forward problem with surface
integrals on the boundaries of the head compartments, meaning that the linear
systems to solve are considerably smaller. The BEM is numerically stable and
can be further augmented using acceleration techniques such as the Fast
Multipole Method (FMM) to reach great computational efficiency [22, 23, 24].
Unfortunately, the standard BEM requires the head medium to be piecewise-
uniform, and is consequently unable to model the anisotropies and
inhomogeneities of the skull and of the white matter. Therefore, this
intrinsic assumption drastically reduces its applicability.

A BEM formulation that accounts for the skull anisotropy has recently been
proposed [25] but relies on a volumetric discretization of the whole head. On
the other hand, the white matter fibers can be quantitatively reconstructed
from tractography algorithms applied on dMRI data [26, 27]. This
representation is promising in that it describes white matter in a much more
effective and precise way than when relying on a complete volume
discretization which is independent of the fiber structure. The thin fiber
geometry is a common structure in high frequency electromagnetic problems [28]
and has been studied in a few contributions for bioelectromagnetic problems
[29, 30, 31, 32]. Tractography has also been used in transcranial brain
stimulation studies [33, 34], although the fiber structure was neglected in
the forward model construction. At present and to the best of our knowledge,
no contribution has been proposed to solve the anisotropic forward problem
without resorting to a full volumetric discretization of the head.

In this paper, we address this issue by modifying and complementing the
standard BEM (2D) equations with adequate anisotropy-handling wire (1D)
equations for the white matter and volume (3D) equations for the skull. The
resulting new hybrid integral formulation is effectively tailored to the EEG
forward problem as every head tissue is suitably discretized according to its
electrical properties. The integration of tractography algorithms for white
matter conductivity profiling is obtained with an electrically coherent
derivation of the fiber parameters, enabling the computation of a multimodal
MRI consistent solution of the anisotropic EEG forward problem. The validity
of the proposed scheme is confirmed by numerical experiments which demonstrate
its practical relevance. Very preliminary results from this work were
presented at the conferences [35, 36].

The background and notation required throughout this paper are introduced in
Section 2 and the new formulation is developed along with its discretization
in Section 3. Diffusion MRI considerations for white matter modeling are
discussed in Section 4, followed by numerical results in Section 5. Conclusive
remarks are presented in Section 6.

##  2 Background and Notations

In the following, the head medium $\Omega$ will be modeled by nested layers
$\Omega_{i}$ ($\Omega=\bigcup_{i=1}^{N}\Omega_{i}$) with respective external
boundaries $\Gamma_{i}$ ($\Gamma=\bigcup_{i=1}^{N}\Gamma_{i})$ and scalar
background conductivities $\sigma_{i}$, typically accounting for brain, skull
and scalp [21]. We assume a nested structure for the sake of simplicity and
readability but the formulation can be extended to arbitrary geometries. In
all generality, the conductivity is expressed as a symmetric positive
semidefinite spatial tensor

$\bar{\bar{\sigma}}(\bm{r})=U(\bm{r})\begin{bmatrix}\sigma_{v_{1}}(\bm{r})&0&0\\\
0&\sigma_{v_{2}}(\bm{r})&0\\\
0&0&\sigma_{v_{3}}(\bm{r})\end{bmatrix}U^{T}(\bm{r}),$ (1)

where $\sigma_{v_{k}}$ are the conductivity eigenvalues and the unitary matrix
$U$ contains the associated eigenvectors. Classical BEM assumes that in each
compartment, the conductivity is reasonably approximated by the background
conductivity so that
$\bar{\bar{\sigma}}(\bm{r}\in\Omega_{i})\approx\sigma_{i}\mathbfsf{I}$, where
$\mathbfsf{I}$ is the identity matrix. In practice, this is poorly verified in
the skull, and in the fibrous white matter compartment $\Omega_{i_{w}}$, in
which ions flow preferentially along the local fiber direction
$\hat{\bm{l}}(\bm{r})$ so that
$\bar{\bar{\sigma}}=\sigma_{i_{w}}\mathbfsf{I}+(\sigma_{l}-\sigma_{i_{w}})\hat{\bm{l}}\hat{\bm{l}}^{T}\,,$
where the longitudinal conductivity $\sigma_{l}$ is approximately ten times
higher than the background conductivity $\sigma_{i_{w}}$ [37, 38]. The
external region $\Omega_{N+1}$ represents air and it is not conductive,
therefore $\sigma_{N+1}=0$. In this framework, solving the forward problem
amounts to determining the unknown electric potential $\phi(\bm{r})$ on the
scalp surface $\Gamma_{N}$ generated by a known source located at
$\bm{r}_{0}\in\Omega$ inducing a primary current $\bm{J}_{p}(\bm{r})$. In the
quasi-static regime, these quantities are linked by Poisson’s equation

$\nabla\cdot\left(\bar{\bar{\sigma}}\nabla\phi\right)=\nabla\cdot\bm{J}_{s},\quad\bm{r}\in\Omega\,,$
(2)

with the boundary conditions

$\displaystyle\phi|^{-}_{i}$
$\displaystyle=\phi|^{+}_{i}\,,\quad\bm{r}\in\Gamma_{i}\,,$ (3a)
$\displaystyle\hat{\bm{n}}\cdot\bar{\bar{\sigma}}\nabla\phi|^{-}_{i}$
$\displaystyle=\hat{\bm{n}}\cdot\bar{\bar{\sigma}}\nabla\phi|^{+}_{i}\,,\quad\bm{r}\in\Gamma_{i<N}\,,$
(3b) $\displaystyle\hat{\bm{n}}\cdot\bar{\bar{\sigma}}\nabla\phi$
$\displaystyle=0\,,\quad\bm{r}\in\Gamma_{N}\,,$ (3c)

where the dependency in $\bm{r}$ has been dropped for the sake of brevity,
$f|^{+}_{i}$ and $f|^{-}_{i}$ denote respectively the exterior and interior
limits of a function $f$ on surface $\Gamma_{i}$, $\hat{\bm{n}}$ is the
unitary outward normal vector of $\Gamma_{i}$ and $\bm{J}_{s}$ is the primary
source term; in the current setting $\bm{J}_{s}=\bm{J}_{p}$. In the standard
single-layer BEM framework [18], (2) and (3) are reformulated in terms of
integral equations and we define the surface and volume operators

$\displaystyle(\mathcal{S}f_{s})(\bm{r})$
$\displaystyle=\int_{S}G\left(\bm{r},\bm{r}^{\prime}\right)f_{s}\left(\bm{r}^{\prime}\right)\mathrm{d}S^{\prime}\,,\quad\bm{r}\in\Omega\,,$
(4a) $\displaystyle(\mathcal{D}^{*}f_{s})(\bm{r})$
$\displaystyle=\int_{S}\nabla_{n}G(\bm{r},\bm{r}^{\prime})f_{s}(\bm{r})\mathrm{d}S^{\prime}\,,\quad\bm{r}\in\Gamma\,,$
(4b) $\displaystyle(\mathcal{S}^{*}_{v}\bm{f}_{v})(\bm{r})$
$\displaystyle=\int_{V}G(\bm{r},\bm{r}^{\prime})\nabla^{\prime}\cdot\bm{f}_{v}(\bm{r}^{\prime})\mathrm{d}V^{\prime}\,,\quad\bm{r}\in\Omega\,,$
(4c) $\displaystyle(\mathcal{D}^{*}_{v}\bm{f}_{v})(\bm{r})$
$\displaystyle=\int_{V}\nabla_{n}G(\bm{r},\bm{r}^{\prime})\nabla^{\prime}\cdot\bm{f}_{v}\left(\bm{r}^{\prime}\right)\mathrm{d}V^{\prime}\,,\quad\bm{r}\in\Gamma\,,$
(4d)

where $\nabla_{n}=\hat{\bm{n}}\cdot\nabla$, $f_{s}$ (resp. $\bm{f}_{v}$) is a
scalar (resp. vector) function defined on a surface $S$ (resp. volume $V$),
and $G(\bm{r},\bm{r}^{\prime})=1/4\pi\|\bm{r}-\bm{r}^{\prime}\|$ is the static
Green’s function. To describe the disparity between background isotropic and
actual conductivity, we define in each compartment $\Omega_{i}$ the
conductivity contrast
$\bar{\bar{\chi}}_{i}(\bm{r})=(\sigma_{i}\mathbfsf{I}-\bar{\bar{\sigma}}(\bm{r}))\bar{\bar{\sigma}}^{-1}(\bm{r})$.
The exterior domain is uniform, so we can write $\bar{\bar{\chi}}_{N+1}=0$
without loss of generality. Since the conductivity contrast is not negligible
in realistic skull and white matter models, the applicability of a pure
surface integral equation method is compromised.

##  3 A New Hybrid Integral Equation Formulation

###  3.1 Integral equations for surface contributions

In all generality, Poisson’s equation (2) in each compartment can be rewritten

$\sigma_{i}\Delta\phi(\bm{r})=\nabla\cdot(\bm{J}_{p}(\bm{r})+\bar{\bar{\chi}}_{i}(\bm{r})\bar{\bar{\sigma}}(\bm{r})\nabla\phi(\bm{r})),\quad\bm{r}\in\Omega_{i}.$
(5)

Hence, a non-uniform conductivity in $\Omega_{i}$ is equivalent to a
volumetric source generating a current
$\bm{J}_{eq_{i}}=\bar{\bar{\chi}}_{i}\bar{\bar{\sigma}}\nabla\phi$ in a medium
with arbitrary uniform background conductivity $\sigma_{i}$. The left-hand
side of (5) implies that the anisotropic forward problem can be recasted as a
piecewise-homogeneous one (with conductivity $\sigma_{i}$ in each compartment)
provided that the source term $\bm{J}_{s}=\bm{J}_{p}+\sum_{k}\bm{J}_{eq_{k}}$
now contains the equivalent current contributions. Therefore, following the
classical BEM [18], we define the single-layer potential

$\phi_{SL}=\sum_{k=1}^{N}\mathcal{S}\xi_{\Gamma_{k}}=\phi+\tfrac{1}{\sigma_{p}}\mathcal{S}^{*}_{v}\bm{J}_{p}+\sum_{k=1}^{N}\tfrac{1}{\sigma_{k}}\mathcal{S}^{*}_{v}\bm{J}_{eq_{k}}\,,$
(6)

where $\xi_{\Gamma_{i}}=[\nabla_{n}\phi]_{\Gamma_{i}}$. On each interface
$\Gamma_{i}$, the boundary condition (3b) imposes

$\displaystyle[\hat{\bm{n}}\cdot\bar{\bar{\sigma}}\nabla\phi]_{\Gamma_{i}}=[\sigma_{i}\nabla_{n}\phi]_{\Gamma_{i}}-[\hat{\bm{n}}\cdot\bm{J}_{eq_{i}}]_{\Gamma_{i}}=0.$
(7)

Using (6), the term $[\sigma_{i}\nabla_{n}\phi]_{\Gamma_{i}}$ can be broken
into

$\displaystyle[\sigma_{i}\nabla_{n}\phi_{SL}]_{\Gamma_{i}}$
$\displaystyle=\tfrac{\sigma_{i}+\sigma_{i+1}}{2}\xi_{\Gamma_{i}}+(\sigma_{i}-\sigma_{i+1})\sum_{k}\mathcal{D}^{*}\xi_{\Gamma_{k}}\,,$
$\displaystyle[\sigma_{i}\nabla_{n}\tfrac{1}{\sigma_{p}}\mathcal{S}^{*}_{v}\bm{J}_{p}]_{\Gamma_{i}}$
$\displaystyle=\tfrac{\sigma_{i}-\sigma_{i+1}}{\sigma_{p}}\mathcal{D}^{*}_{v}\bm{J}_{p}\,,$
$\displaystyle\sum_{k}[\tfrac{\sigma_{i}}{\sigma_{k}}\nabla_{n}\mathcal{S}^{*}_{v}\bm{J}_{eq_{k}}]_{\Gamma_{i}}$
$\displaystyle=(\sigma_{i}-\sigma_{i+1})\sum_{k}\tfrac{1}{\sigma_{k}}\mathcal{{D}}^{*}_{v}\bm{J}_{eq_{k}}$
$\displaystyle+\hat{\bm{n}}\cdot$
$\displaystyle(\tfrac{\sigma_{i}+\sigma_{i+1}}{2\sigma_{i}}\bm{J}_{eq_{i}}-\tfrac{\sigma_{i}+\sigma_{i+1}}{2\sigma_{i+1}}\bm{J}_{eq_{i+1}})\,,$
(8)

where we used the fact that the conductivity contrast, and therefore the
equivalent currents $\bm{J}_{eq_{i}}$ and $\bm{J}_{eq_{i+1}}$ are
discontinuous across $\Gamma_{i}$, so that the double-layer volume operator
$\mathcal{D}^{*}_{v}$ in (3.1) gives rise to the discontinuity

$\displaystyle[\mathcal{D}^{*}_{v}\bm{J}_{eq_{i}}]_{\Gamma_{i}}$
$\displaystyle=\hat{\bm{n}}\cdot\bm{J}_{eq_{i}}\,,\quad\bm{r}\in\Gamma_{i}\,,$
(9) $\displaystyle[\mathcal{D}^{*}_{v}\bm{J}_{eq_{i+1}}]_{\Gamma_{i}}$
$\displaystyle=-\hat{\bm{n}}\cdot\bm{J}_{eq_{i+1}}\,,\quad\bm{r}\in\Gamma_{i}\,.$
(10)

Substituting (3.1) in (7), the equation on each surface $\Gamma_{i}$ is

$\tfrac{\sigma_{i}+\sigma_{i+1}}{2(\sigma_{i+1}-\sigma_{i})}\xi_{\Gamma_{i}}+\hat{\bm{n}}\cdot(\tfrac{1}{2\sigma_{i+1}}\bm{J}_{eq_{i+1}}-\tfrac{1}{2\sigma_{i}}\bm{J}_{eq_{i}})\\\
-\sum_{k}[\mathcal{D}^{*}\xi_{\Gamma_{k}}-\tfrac{1}{\sigma_{k}}\mathcal{{D}}^{*}_{v}\bm{J}_{eq_{k}}]=-\tfrac{1}{\sigma_{p}}\mathcal{D}^{*}_{v}\bm{J}_{p}\,.$
(11)

To obtain a discrete system, the surface unknowns $\xi_{\Gamma_{i}}$ are
expanded with pyramidal basis functions
$\xi_{\Gamma}(\bm{r})=\sum_{j}\alpha_{p_{j}}p_{j}(\bm{r})$, where $p_{j}$ is
non-zero in any triangle $T_{jkl}$ defined by vertices
$\bm{r}_{j},\bm{r}_{k},\bm{r}_{l}$

$p_{j}(\bm{r})=\begin{cases}\frac{|(\bm{r}-\bm{r}_{l})\times(\bm{r}_{k}-\bm{r}_{l})|}{|(\bm{r}_{j}-\bm{r}_{l})\times(\bm{r}_{k}-\bm{r}_{l})|},&\bm{r}\in
T_{jkl}\,,\\\ 0,&\text{otherwise}.\end{cases}$ (12)

Following a Galerkin approach the surface equations are tested with the same
basis functions, resulting in the adjoint double-layer and sparse Gram self-
terms

$\displaystyle(\mathsfbfit{G}_{pp})_{mn}$
$\displaystyle=\tfrac{\sigma_{i_{n}}+\sigma_{i_{n}+1}}{2(\sigma_{i_{n}+1}-\sigma_{i_{n}})}\langle
p_{m},p_{n}\rangle_{\Gamma},$ (13a)
$\displaystyle(\mathsfbfit{D}^{*}_{pp})_{mn}$ $\displaystyle=\langle
p_{m},\mathcal{D}^{*}p_{n}\rangle_{\Gamma}\,,$ (13b)

in which $\langle
a,b\rangle_{\kappa}=\int_{\kappa}a(\bm{r})b(\bm{r})\mathrm{d}\kappa$, and $m$
and $n$ denote a testing row and a basis column, respectively.

Contrary to $\bm{J}_{p}$, the equivalent currents are secondary sources in
that they arise from the primary activity and are dependent on the media.
Therefore they must be determined by coupling (11) with additional equations.

###  3.2 Integral equations for volume contributions

By applying the gradient operator to (6) in the inhomogeneous regions, we
obtain a new set of independent equations

$-(\sigma_{i}\mathbfsf{I}-\bar{\bar{\sigma}})^{-1}\bm{J}_{eq_{i}}+\sum_{k}\left[\nabla\mathcal{S}\xi_{\Gamma_{k}}-\frac{1}{\sigma_{k}}\nabla\mathcal{S}^{*}_{v}\bm{J}_{eq_{k}}\right]\\\
=\frac{1}{\sigma_{p}}\nabla\mathcal{S}^{*}_{v}\bm{J}_{p}\,.$ (14)

Note that (14) is valid when the conductivity contrast is not null, so that
the volume equations requires the discretization of the inhomogeneous domain
only. The equivalent volume currents are expanded with Schaubert-Wilton-
Glisson (SWG) basis functions [39]
$\bm{J}_{eq}(\bm{r})=\bar{\bar{\chi}}(\bm{r})\sum_{j}\alpha_{s_{j}}\bm{s}_{j}(\bm{r})$,
where $\bm{s}_{j}$ is defined on the tetrahedron pair $T_{j}^{\pm}$ with
volumes $V_{j}^{\pm}$, opposite vertices $\bm{r}_{j}^{\pm}$ and common
triangle $t_{j}$ and area $a_{j}$, as

$\bm{s}_{j}(\bm{r})=\begin{cases}\frac{a_{j}}{3V_{j}^{+}}(\bm{r}-\bm{r}_{j}^{+})\,,&\bm{r}\in
T_{j}^{+},\\\ -\frac{a_{j}}{3V_{j}^{-}}(\bm{r}-\bm{r}_{j}^{-})\,,&\bm{r}\in
T_{j}^{-},\\\ \bm{0}\,,&\text{otherwise}.\end{cases}$ (15)

This choice of basis function enforces the continuity of the current density
across each pair of tetrahedra, and also allows integration by parts to
transfer the source gradient to the testing side when computing the entries of
the volume matrix. The volume self terms are

$\displaystyle(\mathsfbfit{G_{v}}_{ss})_{mn}$
$\displaystyle=\langle\bm{s}_{m},\bar{\bar{\chi}}^{-1}\bm{s}_{n}\rangle_{\Omega},$
(16a) $\displaystyle(\mathsfbfit{S^{*}_{v}}_{ss})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle\bm{s}_{m},\nabla\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\Omega}$
$\displaystyle=-\tfrac{1}{\sigma_{i_{n}}}\langle\nabla\cdot\bm{s}_{m},\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\Omega}\,,$
(16b)

where, for implementation purposes, the differential operator on the source
integral is transferred to the testing side via integration by parts and Gauss
divergence theorem. In particular, on the boundary of the tetrahedral mesh,
the divergence of the half SWG (in both test and source integral) gives rise
to a surface integral on its defining triangle $t_{m}$, e.g.

$\displaystyle({\mathsfbfit{S^{*}_{v}}}_{ss})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle\hat{\bm{n}}\cdot\bm{s}_{m},\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\partial\Omega}\
-\tfrac{1}{\sigma_{i_{n}}}\langle\nabla\cdot\bm{s}_{m},\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\Omega}\
$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\int_{t_{m}}\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}(\bm{r})\mathrm{d}S-\tfrac{a_{m}}{\sigma_{i_{n}}V_{m}^{+}}\int_{T_{m}^{+}}\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}(\bm{r})\mathrm{d}V$
(17)

for a half testing SWG function $\bm{s}_{m}$.

###  3.3 Integral equations for fiber contributions

The white matter compartment $\Omega_{i_{w}}$ is locally more conductive along
a specific direction $\hat{\bm{l}}(\bm{r})$, so that the non-uniform region is
locally modeled as a thin cylindrical fiber. The fiber contrast
$\bar{\bar{\chi}}_{i_{w}}$ is a projection on $\hat{\bm{l}}$, which implies
that the volume unknown along a fiber,
$\bm{J}_{eq_{i_{w}}}=(\sigma_{i_{w}}-\sigma_{l})\frac{\partial\phi}{\partial
l}\hat{\bm{l}}=J_{eq_{i_{w}}}(l)\hat{\bm{l}}\,,$ is a scalar function along
the fiber direction. Hence, by taking the derivative of (6) along
$\hat{\bm{l}}$ and with further manipulations, we obtain the wire integral
equations

$-\frac{1}{\sigma_{i_{w}}-\sigma_{l}}J_{eq_{i_{w}}}+\sum_{k}[\nabla_{l}\mathcal{S}\xi_{\Gamma_{k}}-\tfrac{1}{\sigma_{k}}\nabla_{l}\mathcal{S}^{*}_{v}\bm{J}_{eq_{k}}]\\\
=\tfrac{1}{\sigma_{p}}\nabla_{l}\mathcal{S}^{*}_{v}\bm{J}_{p},\quad\bar{\bar{\chi}}_{i_{w}}(\bm{r})\neq
0\,.$ (18)

Following the fiber conductivity model, the fiber equivalent currents in the
white matter are expanded using one-dimensional basis functions with
longitudinal orientation and constant value on the fiber cross-section.
Piecewise linear (hat) basis functions defined on consecutive fiber segments
$s^{-}=[\bm{r}_{j-1};\bm{r}_{j}]$ and $s^{+}=[\bm{r}_{j};\bm{r}_{j+1}]$ as

$\bm{h}_{j}(\bm{r})=\begin{cases}\frac{1}{|\bm{r}_{j}-\bm{r}_{j-1}|}(\bm{r}-\bm{r}_{j-1})\,,&\bm{r}\in
s^{-},\\\
\frac{1}{|\bm{r}_{j+1}-\bm{r}_{j}|}(\bm{r}_{j+1}-\bm{r})\,,&\bm{r}\in
s^{+},\\\ \bm{0}\,,&\text{otherwise},\end{cases}$ (19)

automatically enforce continuity of the current density across segments and,
similarly to the SWG functions in the 3D case, are divergence conforming,
which allows the transfer of the gradient in the source integral to the
testing function. The fiber self terms are

$\displaystyle(\mathsfbfit{G_{v}}_{hh})_{mn}$
$\displaystyle=\langle\bm{h}_{m},\bar{\bar{\chi}}^{-1}\bm{h}_{n}\rangle_{\Omega},$
(20a) $\displaystyle(\mathsfbfit{S^{*}_{v}}_{hh})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle\bm{h}_{m},\nabla\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{h}_{n}\rangle_{\Omega}\,.$
(20b)

The $\mathsfbfit{S^{*}_{v}}_{hh}$ operator can be efficiently computed with
one-dimensional integrals along the fiber direction or semi-analytically for
far and near interactions, respectively.

###  3.4 Coupling terms

As a direct result of the coupling between the surface, volume and fiber
unknowns in the associated equations (11), (14) and (18), the hybrid system
matrix contains non-zero off-diagonal terms.

$\displaystyle(\mathsfbfit{G}_{ps})_{mn}$
$\displaystyle=\tfrac{1}{2\sigma_{i_{n}}}\langle
p_{m},\hat{\bm{n}}\cdot\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\Gamma},$ (21a)
$\displaystyle(\mathsfbfit{D^{*}_{v}}_{ps})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle
p_{m},\mathcal{D}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\Gamma},$ (21b)
$\displaystyle(\mathsfbfit{D^{*}_{v}}_{ph})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle
p_{m},\mathcal{D}_{v}^{*}\bar{\bar{\chi}}\bm{h}_{n}\rangle_{\Gamma},$ (21c)
$\displaystyle(\mathsfbfit{S}_{sp})_{mn}$
$\displaystyle=\langle\bm{s}_{m},\nabla\mathcal{S}p_{n}\rangle_{\Omega},$
(21d) $\displaystyle(\mathsfbfit{S}_{hp})_{mn}$
$\displaystyle=\langle\bm{h}_{m},\nabla\mathcal{S}p_{n}\rangle_{\Omega},$
(21e) $\displaystyle(\mathsfbfit{S^{*}_{v}}_{hs})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle\bm{h}_{m},\nabla\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{s}_{n}\rangle_{\Omega},$
(21f) $\displaystyle(\mathsfbfit{S^{*}_{v}}_{sh})_{mn}$
$\displaystyle=\tfrac{1}{\sigma_{i_{n}}}\langle\bm{s}_{m},\nabla\mathcal{S}_{v}^{*}\bar{\bar{\chi}}\bm{h}_{n}\rangle_{\Omega}.$
(21g)

The mixed Gram operators between the surface and volume unknowns are a
consequence of the fact that a volume-discretized compartment $\Omega_{i}$ is
bounded by the surface-discretized interfaces $\Gamma_{i}$ and $\Gamma_{i-1}$.
Like the self terms in (16b) and (3.2), the source differential operators are
transferred to the divergence-conforming testing side and computed using
similar strategies as in [40, 41].

###  3.5 Solution of the resulting anisotropic forward problem

The right-hand side (RHS) of the discretized equations is

$\displaystyle(\mathsfbfit{v}_{p})_{m}$
$\displaystyle=-\tfrac{1}{\sigma_{s}}\langle
p_{m},\mathcal{D}^{*}_{v}\bm{J}_{p}\rangle_{\Gamma},$ (22a)
$\displaystyle(\mathsfbfit{v}_{h})_{m}$
$\displaystyle=-\tfrac{1}{\sigma_{p}}\langle\bm{h}_{m},\nabla\mathcal{S}^{*}_{v}\bm{J}_{p}\rangle_{\Omega}=\tfrac{1}{\sigma_{p}}\langle\nabla\cdot\bm{h}_{m},\mathcal{S}^{*}_{v}\bm{J}_{p}\rangle_{\Omega},$
(22b) $\displaystyle(\mathsfbfit{v}_{s})_{m}$
$\displaystyle=-\tfrac{1}{\sigma_{p}}\langle\bm{s}_{m},\nabla\mathcal{S}^{*}_{v}\bm{J}_{p}\rangle_{\Omega}=\tfrac{1}{\sigma_{p}}\langle\nabla\cdot\bm{s}_{m},\mathcal{S}^{*}_{v}\bm{J}_{p}\rangle_{\Omega}.$
(22c)

The hybrid matrix system

$\displaystyle\small{\begin{bmatrix}-\mathsfbfit{G}_{pp}+\mathsfbfit{D}^{*}_{pp}&-\mathsfbfit{D^{*}_{v}}_{ph}&\mathsfbfit{G}_{ps}-\mathsfbfit{D^{*}_{v}}_{ps}\\\
-\mathsfbfit{S}_{hp}&\mathsfbfit{G_{v}}_{hh}+\mathsfbfit{S^{*}_{v}}_{hh}&\mathsfbfit{S^{*}_{v}}_{hs}\\\
-\mathsfbfit{S}_{sp}&\mathsfbfit{S^{*}_{v}}_{sh}&\mathsfbfit{G_{v}}_{ss}+\mathsfbfit{S^{*}_{v}}_{ss}\end{bmatrix}\begin{bmatrix}\alpha_{p}\\\
\alpha_{h}\\\ \alpha_{s}\end{bmatrix}=\begin{bmatrix}\mathsfbfit{v}_{p}\\\
\mathsfbfit{v}_{h}\\\ \mathsfbfit{v}_{s}\end{bmatrix}}$ (23)

has a well-known singularity in the surface-surface (top-left) block which is
removed by deflation [42]. After inverting the final system, the potential in
$\Omega$ is finally computed as

$\phi(\bm{r})=\sum_{i}\alpha_{p_{i}}\mathcal{S}p_{i}(\bm{r})-\sum_{j}\alpha_{h_{j}}\mathcal{S}^{*}_{v}\bar{\bar{\chi}}\bm{h}_{j}(\bm{r})\\\
-\sum_{k}\alpha_{s_{k}}\mathcal{S}^{*}_{v}\bar{\bar{\chi}}\bm{s}_{k}(\bm{r})-\tfrac{1}{\sigma_{p}}\mathcal{S}^{*}_{v}\bm{J}_{p}(\bm{r}).$
(24)

##  4 A New Diffusion MRI Consistent White Matter Model

Although tractography provides the mesh for the fibers, their radius,
longitudinal and background conductivity must be determined. To this end, we
propose next a new procedure to enable and validate the proposed tractography-
derived conductivity model of white matter. First the fiber radius is adjusted
to match a known volume of white matter (approximately
$450\text{\,}{\mathrm{cm}}^{3}$ [43]). Then, for a given fiber conductivity,
an effective homogenized conductivity of the brain $\bar{\bar{\sigma}}_{e}$ is
determined by solving a homogenization problem with integral equations.
Denoting $\Omega_{1}$ the brain compartment (which contains the fibers),
$\bar{\bar{\sigma}}_{e}$ is defined as
$\langle\bm{J}\rangle=\bar{\bar{\sigma}}_{e}\langle\bm{E}\rangle$, where
$\bm{J}$ is the current density, $\bm{E}=-\nabla\phi$ is the electric field
and $\langle\cdot\rangle=V^{-1}_{\Omega_{1}}\int_{\Omega_{1}}{\cdot\
d\Omega_{1}}$ denotes the volume average over the brain. By imposing a normal
current density $\bm{J}_{0}$ on the boundary, that is

$\hat{\bm{n}}\cdot\bar{\bar{\sigma}}\nabla\phi=\hat{\bm{n}}\cdot\bm{J}_{0},\
\bm{r}\in\Gamma_{1},$ (25)

and using Gauss divergence theorem, the volume averaged current density
$\langle\bm{J}\rangle$ is equal to $\bm{J}_{0}$. Furthermore, a forward
problem for the Poisson’s equation in the absence of primary source in the
volume, $\nabla\cdot\bar{\bar{\sigma}}\nabla\phi=0,\ \bm{r}\in\Omega_{1},$ and
complemented with the boundary condition (25) can be solved with our hybrid
formulation to obtain the scalar potential $\phi$ on the boundary
$\Gamma_{1}$, as given by (24). Using again Gauss divergence theorem, the
computation of $\langle\bm{E}\rangle$ is given by

$\langle\bm{E}^{(\bm{J}_{0})}\rangle=-\frac{1}{V_{\Omega_{1}}}\int_{\Gamma_{1}}{\phi(\bm{r})\hat{\bm{n}}(\bm{r})d\Gamma_{1}},$
(26)

where the $\bm{J}_{0}$ superscript indicates the dependency on the imposed
current. Thus, by solving this forward problem three times for
$\bm{J}_{0}\in\\{\bm{x},\bm{y},\bm{z}\\}$ (the unit cartesian vectors), the
equivalent conductivity is obtained as

$\bar{\bar{\sigma}}_{e}=\begin{bmatrix}\langle\bm{E}^{(\bm{x})}\rangle_{x}&0&0\\\
0&\langle\bm{E}^{(\bm{y})}\rangle_{y}&0\\\
0&0&\langle\bm{E}^{(\bm{z})}\rangle_{z}\end{bmatrix}^{-1},$ (27)

where the subscripts on vectors $\langle\bm{E}^{(\bm{J}_{0})}\rangle$ denote
their cartesian components. As $\bar{\bar{\sigma}}_{e}$ satisfies Ohm’s law
averaged over the volume, it corresponds to the electrically equivalent
average (homogenized) conductivity of the brain for the specified boundary
conditions (25). This value can then be compared with brain conductivity
values reported in literature to ensure that the obtained tractography-based
conductivity model of the white matter is realistic. This will be demonstrated
in the numerical results.

(a)

$40$$60$$80$$100$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$Dipole Eccentricity
$\left[\mathrm{$\%$}\right]$Relative Error $\left[\mathrm{$\%$}\right]$BEM
(isotropic)FEMThis work (b)

Figure 1: (a) Spherical head geometry with anisotropic and inhomogeneous brain and skull compartments; (b) relative error of the proposed method, OpenMEEG BEM and FieldTrip FEM. Table 1: Partial and total timing of the proposed method and FEM in a realistic scenario | Setup time | Time per RHS | Total time (4 elec.) | Total time (16 elec.) | Total time (64 elec.) | Total time (256 elec.)
---|---|---|---|---|---|---
This work | $2262\text{\,}\mathrm{s}$ | $0.6\text{\,}\mathrm{s}$ | $2264\text{\,}\mathrm{s}$ | $2271\text{\,}\mathrm{s}$ | $2301\text{\,}\mathrm{s}$ | $2416\text{\,}\mathrm{s}$
FEM | $81\text{\,}\mathrm{s}$ | $74\text{\,}\mathrm{s}$ | $303\text{\,}\mathrm{s}$ | $1191\text{\,}\mathrm{s}$ | $4742\text{\,}\mathrm{s}$ | $18\,950\text{\,}\mathrm{s}$

##  5 Simulation Results

###  5.1 Validation on a canonical model

Following a well established practice in literature, we benchmark the accuracy
of our new formulation with a standard 3-layer (brain, skull, scalp) spherical
model, with normalized radii of $0.87$, $0.92$, $1$ and normalized background
conductivities of $11\text{/}301\text{\,}\mathrm{,}$ respectively. The skull
inhomogeneous conductivity is modeled using the same anisotropic model and as
in [44], with a varying soft bone thickness that is equal to 1 in an arbitrary
position $\bm{p}=[0.9,0,0]^{T}$ and decreases linearly to 0 according to the
distance to $\bm{p}$. The white matter anisotropy is modeled in the first
sphere by adding $30$ fibers with a radius of 0.05, a length of 1.4 and a
longitudinal conductivity that is ten times that of the background. A current
dipole acting as the primary source is placed along the $z$ axis with a dipole
moment of $[1,1,1]$ in the $xyz$-coordinate system. The geometry, discretized
with an average edge length $h=0.1$, is illustrated in Figure 1a. Leveraging
on the reciprocity principle, the hybrid system is directly inverted via LU
factorization and backward substitution. A higher resolution FEM solution with
quadratic basis functions ($3\,110\,264$ unknowns, $h=0.03$) is used as
reference and the accuracy of the hybrid scheme ($23\,358$ unknowns) is
plotted in Figure 1b as the relative $\ell_{2}$-error of the scalp electric
potential, for different dipole excentricities. For comparison, a piecewise-
isotropic symmetric BEM solution ($12\,648$ unknowns) is computed with
OpenMEEG [45], along with an anisotropic FEM solution ($394\,807$ unknowns)
obtained with the FieldTrip toolbox [46]. The comparison with the symmetric
BEM shows how standard surface modeling, which omits our fiber and volume
integral equations, results in important errors (around
$20\text{\,}\mathrm{\char 37\relax}$). In contrast, performing slightly better
than the completely volumetric FEM solution, our 3-2-1D hybrid formulation
exhibits less than $1\text{\,}\mathrm{\char 37\relax}$ relative error up to
very high source excentricity. This experiment confirms the ability of our
hybrid formulation to overcome the model approximations of a standard BEM
approach.

Figure 2: Cortex map of the relative error between the anisotropic FEM and
hybrid solver. Figure 3: (a) Skull volume currents, (b) scalp surface
potential and (c) fiber currents computed with the hybrid formulation. The
arrow in (c) represents the cortical dipole and the black dots in (b) the
electrode positions.

###  5.2 Validation of the conductivity model

The hybrid solver is also tested and compared with an anisotropic FEM solution
on a realistic head model derived from the MRI of a single subject from the
Wu-Minn Human Connectome Project database [47] to highlight its applicability.
Surface and volume meshes were obtained after preprocessing [48], segmentation
and tessellation [49] of the structural MRI data. The different head
compartments were assigned commonly used conductivity values [38]. The skull
anisotropy is modeled following the model of [44], with soft bone thickness
defined as an affine function of the local skull thickness. A probabilistic
tractography algorithm [50] was applied on the subject’s diffusion MRI data to
generate a bundle of non-connected streamlines, which was subsequently
clustered [51] to obtain a concise fiber map of the white matter. The fiber
radius was adjusted to match a white matter volume of
$450\text{\,}{\mathrm{cm}}^{3}$. We applied the procedure described in Section
4 with the anisotropic longitudinal and transversal conductivities reported in
[37] and obtained the homogenized brain conductivity

$\bar{\bar{\sigma}}_{e}=\begin{bmatrix}0.1755&0&0\\\ 0&0.1915&0\\\
0&0&0.1855\end{bmatrix}\,.$ (28)

These values are fairly close to the isotropic brain conductivity
$\sigma_{brain}=$ $0.18$ $\text{\,}\mathrm{S}\text{\,}{\mathrm{m}}^{-1}$ [38],
and thus confirms that our fiber conductivity model is consistent with the
conductivity values reported in literature.

$0$$20$$40$$60$$80$$-200$$0$Potential
$\left[\mathrm{$\mathrm{\SIUnitSymbolMicro V}$}\right]$Hybrid iso.Hybrid
aniso.FEM iso.FEM aniso.

(a)

$0$$20$$40$$60$$80$$0$$200$$400$$600$Electrode index Potential
$\left[\mathrm{$\mathrm{\SIUnitSymbolMicro V}$}\right]$Hybrid iso.Hybrid
aniso.FEM iso.FEM aniso. (b)

Figure 4: Scalp electrode potential resulting from (a) a deep and (b) a
shallow source obtained from solving the EEG forward problem with the hybrid
integral method (continuous lines) and FEM (dots) and with isotropic (blue)
and anisotropic (red) conductivity assumptions.

###  5.3 Comparison with an anisotropic full volume FEM model

An anisotropic FEM solution was obtained from the same dMRI data by computing
the diffusion tensor via a least-square fit for each white matter voxel. The
conductivity tensor was then derived by following a mean conductivity volume
constraint [52]. Figure 2 displays the difference between our hybrid and the
FEM model, computed as the relative error of the electric potential obtained
on a standardized set of $76$ electrodes (depicted in Figure 3), for each
source on a cortical surface. The forward solutions produced by both numerical
methods are overall in good agreement, although sources in the cortex sulci
exhibit more differences than the shallow gyral sources. This is likely due to
the fact that deep sources are more sensitive to the white matter anisotropy,
which is modeled differently for both methods. This is further illustrated in
Figure 4 which depicts the scalp electrode potential for a deep and a shallow
brain source. The FEM and hybrid forward solution show similar topological and
magnitude changes when accounting for the anisotropy. This numerical
experiment confirms that the proposed multimodal MRI-based hybrid integral
method is consistent with a DTI volume-based anisotropic model of white matter
despite their intrinsic modeling differences.

###  5.4 Timing comparisons

All simulations were run on a 24-core 3.00 GHz Intel Xeon E5-2687W v4 CPU with
768 GB of RAM. Table 1 reports the run times of the two methods for different
electrode configurations. In total, the FEM system contained $1\,025\,411$
unknowns whereas the hybrid solver had $46\,695$. The proposed method incurs a
bigger computational overhead as it must fill an $N^{2}$ (where $N$ is the
number of unknowns) full matrix as opposed to the sparse FEM. Using the
reciprocity principle [53], for $N_{e}$ electrodes, the computation of the
forward matrix requires the solution of $N_{e}-1$ RHS. The computational gain
becomes apparent for medium to high resolution sensor configurations. Like
classical BEM, the $\mathcal{O}(N^{2})$ asymptotic complexity in matrix
building and storage of the proposed integral equation-based method can
furthermore be reduced to linear complexity with fast solvers [54, 22].

##  6 Conclusion

In this paper we have presented a new solution to the anisotropic EEG forward
problem that does not require a full volumetric discretization of the head.
The standard boundary integral formulation was coupled with thin volume and
wire integral equations that adequately match the non-uniform conductivities
of the skull and white matter, respectively. The accuracy of the new BEM-like
and anisotropy-adapted formulation was demonstrated on a canonical model while
a realistic scenario illustrated its applicability in a clinical environment,
in which the patient-specific physiological properties were derived from
advanced biomedical imaging techniques.

## References

  * [1] G. Pfurtscheller, C. Brunner, A. Schlögl, and F. L. Da Silva, “Mu rhythm (de) synchronization and eeg single-trial classification of different motor imagery tasks,” NeuroImage, vol. 31, no. 1, pp. 153–159, 2006.
  * [2] S. Smith, “Eeg in the diagnosis, classification, and management of patients with epilepsy,” Journal of Neurology, Neurosurgery & Psychiatry, vol. 76, no. suppl 2, pp. ii2–ii7, 2005.
  * [3] K. Wendel, O. Väisänen, J. Malmivuo, N. G. Gencer, B. Vanrumste, P. Durka, R. Magjarevic, S. Supek, M. L. Pascu, H. Fontenelle, et al., “Eeg/meg source imaging: methods, challenges, and open issues,” Computational intelligence and neuroscience, vol. 2009, p. 13, 2009.
  * [4] R. D. Pascual-Marqui, “Review of methods for solving the eeg inverse problem,” International journal of bioelectromagnetism, vol. 1, no. 1, pp. 75–86, 1999.
  * [5] S. J. Kiebel, O. David, and K. J. Friston, “Dynamic causal modelling of evoked responses in eeg/meg with lead field parameterization,” NeuroImage, vol. 30, no. 4, pp. 1273–1284, 2006.
  * [6] R. Grech, T. Cassar, J. Muscat, K. P. Camilleri, S. G. Fabri, M. Zervakis, P. Xanthopoulos, V. Sakkalis, and B. Vanrumste, “Review on solving the inverse problem in eeg source analysis,” Journal of neuroengineering and rehabilitation, vol. 5, no. 1, p. 25, 2008.
  * [7] R. Sadleir and A. Argibay, “Modeling skull electrical properties,” Annals of Biomedical engineering, vol. 35, no. 10, pp. 1699–1712, 2007.
  * [8] C. Papageorgakis, Patient specific conductivity models: characterization of the skull bones. PhD thesis, Universite Cote d Azur, 2017.
  * [9] D. Le Bihan and H. Johansen-Berg, “Diffusion mri at 25: exploring brain tissue structure and function,” Neuroimage, vol. 61, no. 2, pp. 324–341, 2012\.
  * [10] W. H. Lee, Z.-D. Deng, T.-S. Kim, A. F. Laine, S. H. Lisanby, and A. V. Peterchev, “Regional electric field induced by electroconvulsive therapy in a realistic finite element head model: influence of white matter anisotropic conductivity,” Neuroimage, vol. 59, no. 3, pp. 2110–2123, 2012.
  * [11] S. Shahid, P. Wen, and T. Ahfock, “Numerical investigation of white matter anisotropic conductivity in defining current distribution under tdcs,” Computer methods and programs in biomedicine, vol. 109, no. 1, pp. 48–64, 2013\.
  * [12] J. Vorwerk, J.-H. Cho, S. Rampp, H. Hamer, T. R. Knösche, and C. H. Wolters, “A guideline for head volume conductor modeling in eeg and meg,” NeuroImage, vol. 100, pp. 590–607, 2014.
  * [13] H. Hallez, B. Vanrumste, R. Grech, J. Muscat, W. De Clercq, A. Vergult, Y. D’Asseler, K. P. Camilleri, S. G. Fabri, S. Van Huffel, et al., “Review on solving the forward problem in eeg source analysis,” Journal of neuroengineering and rehabilitation, vol. 4, no. 1, p. 46, 2007.
  * [14] P. Bruno, J. Hyttinen, P. Inchingolo, A. Magrofuoco, S. Mininel, and F. Vatta, “A fdm anisotropic formulation for eeg simulation,” in Engineering in Medicine and Biology Society, 2006. EMBS’06. 28th Annual International Conference of the IEEE, pp. 1121–1125, IEEE, 2006.
  * [15] V. Montes-Restrepo, P. van Mierlo, G. Strobbe, S. Staelens, S. Vandenberghe, and H. Hallez, “Influence of skull modeling approaches on eeg source localization,” Brain topography, vol. 27, no. 1, pp. 95–111, 2014.
  * [16] S. Lew, C. Wolters, T. Dierkes, C. Röer, and R. MacLeod, “Accuracy and run-time comparison for different potential approaches and iterative solvers in finite element method based eeg source analysis,” Applied Numerical Mathematics, vol. 59, no. 8, pp. 1970–1988, 2009.
  * [17] M. S. Hamalainen and J. Sarvas, “Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data,” IEEE transactions on biomedical engineering, vol. 36, no. 2, pp. 165–171, 1989.
  * [18] J. Kybic, M. Clerc, T. Abboud, O. Faugeras, R. Keriven, and T. Papadopoulo, “A common formalism for the integral formulations of the forward eeg problem,” IEEE transactions on medical imaging, vol. 24, no. 1, pp. 12–28, 2005.
  * [19] L. J. Gomez, A. C. Yücel, and E. Michielssen, “The icvsie: A general purpose integral equation method for bio-electromagnetic analysis,” IEEE Transactions on Biomedical Engineering, vol. 65, no. 3, pp. 565–574, 2017\.
  * [20] E. Niedermeyer and F. L. da Silva, Electroencephalography: basic principles, clinical applications, and related fields.  Lippincott Williams & Wilkins, 2005.
  * [21] R. Brette and A. Destexhe, Handbook of neural activity measurement.  Cambridge University Press, 2012.
  * [22] Y. Liu, Fast multipole boundary element method: theory and applications in engineering.  Cambridge university press, 2009.
  * [23] S. N. Makarov, G. M. Noetscher, T. Raij, and A. Nummenmaa, “A quasi-static boundary element approach with fast multipole acceleration for high-resolution bioelectromagnetic models,” IEEE Transactions on Biomedical Engineering, vol. 65, no. 12, pp. 2675–2683, 2018.
  * [24] A. T. Htet, G. B. Saturnino, E. H. Burnham, G. Noetscher, A. Nummenmaa, and S. N. Makarov, “Comparative performance of the finite element method and the boundary element fast multipole method for problems mimicking transcranial magnetic stimulation (tms),” Journal of neural engineering, 2019.
  * [25] L. Rahmouni, R. Mitharwal, and F. P. Andriulli, “Two volume integral equations for the inhomogeneous and anisotropic forward problem in electroencephalography,” Journal of Computational Physics, vol. 348, pp. 732–743, 2017.
  * [26] K. Yamada, K. Sakai, K. Akazawa, S. I. Yuen, and T. Nishimura, “Mr tractography: a review of its clinical applications.,” Magnetic resonance in medical sciences : MRMS : an official journal of Japan Society of Magnetic Resonance in Medicine, vol. 8 4, pp. 165–74, 2009.
  * [27] B. A. Wandell, “Clarifying human white matter,” Annual review of neuroscience, vol. 39, pp. 103–128, 2016.
  * [28] D. R. Wilton and N. J. Champagne, “Evaluation and integration of the thin wire kernel,” IEEE transactions on antennas and propagation, vol. 54, no. 4, pp. 1200–1206, 2006.
  * [29] E. Olivi, T. Papadopoulo, and M. Clerc, “Handling white-matter anisotropy in bem for the eeg forward problem,” in Proceedings of the 8th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, ISBI 2011., pp. 799–802, IEEE, 2011.
  * [30] A. Pillain, Line, Surface, and Volume Integral Equations for the Electromagnetic Modelling of the Electroencephalography Forward Problem.  PhD thesis, Ecole Nationale Supérieure des Télécommunications de Bretagne-ENSTB, 2016.
  * [31] S. N. Makarov, A. Pascual-Leone, and A. Nummenmaa, “Modeling fiber-like conductivity structures via the boundary element method using thin-wire approximation. i construction of basis functions,” in 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pp. 6473–6476, IEEE, 2016.
  * [32] L. Rahmouni, A. Merlini, A. Pillain, and F. P. Andriulli, “On the boundary element modeling of brain fibers in the eeg forward problem via a new family of wire integral equations,” arXiv preprint arXiv:1903.08414, 2019.
  * [33] A. Nummenmaa, J. A. McNab, P. Savadjiev, Y. Okada, M. S. Hämäläinen, R. Wang, L. L. Wald, A. Pascual-Leone, V. J. Wedeen, and T. Raij, “Targeting of white matter tracts with transcranial magnetic stimulation,” Brain stimulation, vol. 7, no. 1, pp. 80–84, 2014.
  * [34] A. Opitz, M. Windhoff, R. M. Heidemann, R. Turner, and A. Thielscher, “How the brain tissue shapes the electric field induced by transcranial magnetic stimulation,” Neuroimage, vol. 58, no. 3, pp. 849–859, 2011.
  * [35] M. Y. Monin, L. Rahmouni, and F. P. Andriulli, “A hybrid integral equation approach to solve the anisotropic forward problem in electroencephalography,” in 2018 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting, pp. 2395–2396, IEEE, 2018.
  * [36] M. Monin, L. Rahmouni, and F. P. Andriulli, “Diffusion mri consistent wire models for efficient solutions of the anisotropic forward problem in electroencephalography,” in to be presented at 2019 International Conference on Electromagnetics in Advanced Applications (ICEAA), September 2019, 2019.
  * [37] P. W. Nicholson, “Specific impedance of cerebral white matter,” Experimental Neurology, vol. 13, no. 4, pp. 386 – 401, 1965.
  * [38] A. Lozano and M. Hallett, Brain Stimulation.  Handbook of Clinical Neurology, Elsevier Science, 2013.
  * [39] D. Schaubert, D. Wilton, and A. Glisson, “A tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies,” IEEE Transactions on Antennas and Propagation, vol. 32, no. 1, pp. 77–85, 1984.
  * [40] R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-d green’s function or its gradient on a plane triangle,” IEEE Transactions on Antennas and Propagation, vol. 41, no. 10, pp. 1448–1455, 1993.
  * [41] S. Järvenpää, M. Taskinen, and P. Ylä-Oijala, “Singularity extraction technique for integral equation methods with higher order basis functions on plane triangles and tetrahedra,” International journal for numerical methods in engineering, vol. 58, no. 8, pp. 1149–1165, 2003.
  * [42] T. F. Chan, “Deflated decomposition of solutions of nearly singular systems,” SIAM journal on numerical analysis, vol. 21, no. 4, pp. 738–754, 1984.
  * [43] K. Zhang and T. J. Sejnowski, “A universal scaling law between gray matter and white matter of cerebral cortex,” Proceedings of the National Academy of Sciences, vol. 97, no. 10, pp. 5621–5626, 2000.
  * [44] M. Dannhauer, B. Lanfer, C. H. Wolters, and T. R. Knösche, “Modeling of the human skull in eeg source analysis,” Human brain mapping, vol. 32, no. 9, pp. 1383–1399, 2011.
  * [45] A. Gramfort, T. Papadopoulo, E. Olivi, and M. Clerc, “Openmeeg: opensource software for quasistatic bioelectromagnetics,” Biomedical engineering online, vol. 9, no. 1, p. 45, 2010.
  * [46] R. Oostenveld, P. Fries, E. Maris, and J.-M. Schoffelen, “Fieldtrip: open source software for advanced analysis of meg, eeg, and invasive electrophysiological data,” Computational intelligence and neuroscience, vol. 2011, p. 1, 2011.
  * [47] D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. Behrens, E. Yacoub, K. Ugurbil, W.-M. H. Consortium, et al., “The wu-minn human connectome project: an overview,” Neuroimage, vol. 80, pp. 62–79, 2013.
  * [48] M. F. Glasser, S. N. Sotiropoulos, J. A. Wilson, T. S. Coalson, B. Fischl, J. L. Andersson, J. Xu, S. Jbabdi, M. Webster, J. R. Polimeni, et al., “The minimal preprocessing pipelines for the human connectome project,” Neuroimage, vol. 80, pp. 105–124, 2013.
  * [49] A. Thielscher, A. Antunes, and G. B. Saturnino, “Field modeling for transcranial magnetic stimulation: a useful tool to understand the physiological effects of tms?,” in 2015 37th annual international conference of the IEEE engineering in medicine and biology society (EMBC), pp. 222–225, IEEE, 2015.
  * [50] J.-D. Tournier, F. Calamante, and A. Connelly, “Mrtrix: diffusion tractography in crossing fiber regions,” International Journal of Imaging Systems and Technology, vol. 22, no. 1, pp. 53–66, 2012.
  * [51] E. Garyfallidis, M. Brett, M. M. Correia, G. B. Williams, and I. Nimmo-Smith, “Quickbundles, a method for tractography simplification,” Frontiers in neuroscience, vol. 6, p. 175, 2012.
  * [52] M. Rullmann, A. Anwander, M. Dannhauer, S. K. Warfield, F. H. Duffy, and C. H. Wolters, “Eeg source analysis of epileptiform activity using a 1 mm anisotropic hexahedra finite element head model,” NeuroImage, vol. 44, no. 2, pp. 399–410, 2009.
  * [53] C. H. Wolters, L. Grasedyck, and W. Hackbusch, “Efficient computation of lead field bases and influence matrix for the fem-based eeg and meg inverse problem,” Inverse problems, vol. 20, no. 4, p. 1099, 2004.
  * [54] M. Bebendorf, “Approximation of boundary element matrices,” Numerische Mathematik, vol. 86, no. 4, pp. 565–589, 2000.

[]Maxime Y. Monin (S ’18) received the M.Sc. degree in Information
Technologies from IMT Atlantique, Brest, France, and is currently pursuing a
Ph.D. with the Department of Electronics and Telecommunications at the
Politecnico di Torino, Turin, Italy. His research interests are in
computational electromagnetics and their applications for brain imaging.

[]Lyes Rahmouni received the B.Sc. in aeronautical engineering and M.Sc. in
electronic engineering from the University of Science and Technology Houari
Boumediene, Algeria. He has completed a Ph.D. program in computational
electromagnetics from IMT Atlantique, France in 2016, followed by two years of
postdoctoral position at the same grande école. He is currently holder of a
Research fellowship at Politecnico di Torino, Italy. His research interests
focus on integral equation methods and their fast solutions.

[]Adrien Merlini (S ’16, M ’19) received the M.Sc. Eng. degree from the École
Nationale Supérieure des Télécommunications de Bretagne (Télécom Bretagne),
France, in 2015 and received the Ph.D. degree from the École Nationale
Supérieure Mines-Télécom Atlantique (IMT Atlantique), France, in 2019. Between
2018 and 2019, he was a visiting Ph.D. student at the Politecnico di Torino,
Italy, which he then joined as a Research Associate. Since the end of 2019, he
is an Associate Professor with the Microwave Department of IMT Atlantique,
Brest, France.

His research interests include preconditioning and acceleration of integral
equation solvers for electromagnetic simulations and their application in
brain imaging.

[]Francesco P. Andriulli (S ’05, M ’09, SM ’11) received the Laurea in
electrical engineering from the Politecnico di Torino, Italy, in 2004, the MSc
in electrical engineering and computer science from the University of Illinois
at Chicago in 2004, and the PhD in electrical engineering from the University
of Michigan at Ann Arbor in 2008. From 2008 to 2010 he was a Research
Associate with the Politecnico di Torino. From 2010 to 2017 he was an
Associate Professor (2010-2014) and then Full Professor with the École
Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique, previously ENST
Bretagne), Brest, France. Since 2017 he has been a Full Professor with the
Politecnico di Torino, Turin, Italy. His research interests are in
computational electromagnetics with focus on frequency- and time-domain
integral equation solvers, well-conditioned formulations, fast solvers, low-
frequency electromagnetic analyses, and modeling techniques for antennas,
wireless components, microwave circuits, and biomedical applications with a
special focus on Brain Imaging.

Prof. Andriulli was the recipient of the best student paper award at the 2007
URSI North American Radio Science Meeting. He received the first place prize
of the student paper context of the 2008 IEEE Antennas and Propagation Society
International Symposium. He was the recipient of the 2009 RMTG Award for
junior researchers and was awarded two URSI Young Scientist Awards at the
International Symposium on Electromagnetic Theory in 2010 and 2013 where he
was also awarded the second prize in the best paper contest. He also received
the 2015 ICEAA IEEE-APWC Best Paper Award. In addition, he co-authored with
his students and collaborators other three first prize conference papers (EMTS
2016, URSI-DE Meeting 2014, ICEAA 2009), a second prize conference paper (URSI
GASS 2014), a third prize conference paper (IEEE-APS 2018), two honorable
mention conference papers (ICEAA 2011, URSI/IEEE-APS 2013) and other three
finalist conference papers (URSI/IEEE-APS 2012, URSI/IEEE-APS 2007, URSI/IEEE-
APS 2006). Moreover, he received the 2014 IEEE AP-S Donald G. Dudley Jr.
Undergraduate Teaching Award, the triennium 2014-2016 URSI Issac Koga Gold
Medal, and the 2015 L. B. Felsen Award for Excellence in Electrodynamics.

Prof. Andriulli is a member of Eta Kappa Nu, Tau Beta Pi, Phi Kappa Phi, and
of the International Union of Radio Science (URSI). He is the Editor in Chief
of the IEEE Antennas and Propagation Magazine, he serves as a Track Editor for
the IEEE Transactions on Antennas and Propagation, and as an Associate Editor
for the IEEE Antennas and Wireless Propagation Letters, IEEE Access, URSI
Radio Science Letters and IET-MAP. He is the PI of the ERC Consolidator Grant
“321”.