# The Localized Subtraction Approach For EEG and MEG Forward Modeling

Malte B. Höltershinken<EMAIL_ADDRESS>Pia Lange  Tim Erdbrügger
Yvonne Buschermöhle  Fabrice Wallois  Alena Buyx  Sampsa Pursiainen  Johannes
Vorwerk  Christian Engwer  Carsten H. Wolters

###### Abstract

In FEM-based EEG and MEG source analysis, the subtraction approach has been
proposed to simulate sensor measurements generated by neural activity. While
this approach possesses a rigorous foundation and produces accurate results,
its major downside is that it is computationally prohibitively expensive in
practical applications. To overcome this, we developed a new approach, called
the localized subtraction approach. This approach is designed to preserve the
mathematical foundation of the subtraction approach, while also leading to
sparse right-hand sides in the FEM formulation, making it efficiently
computable. We achieve this by introducing a cut-off into the subtraction,
restricting its influence to the immediate neighborhood of the source. In this
work, this approach will be presented, analyzed, and compared to other state-
of-the-art FEM right-hand side approaches. Furthermore, we discuss how to
arrive at an efficient and stable implementation. We perform validation in
multi-layer sphere models where analytical solutions exist. There, we
demonstrate that the localized subtraction approach is vastly more efficient
than the subtraction approach. Moreover, we find that for the EEG forward
problem, the localized subtraction approach is less dependent on the global
structure of the FEM mesh, while the subtraction approach is sensitive to the
mesh structure even at a distance to the source. Additionally, we show the
localized subtraction approach to rival, and in many cases even surpass, the
other investigated approaches in terms of accuracy. Especially for the MEG
forward problem, where Biot-Savart’s law is employed, we show the localized
subtraction approach and the subtraction approach to produce highly accurate
approximations of the volume currents close to the source. The localized
subtraction approach thus reduces the computational cost of the subtraction
approach to an extent that makes it usable in practical applications without
sacrificing the rigorousness and accuracy the subtraction approach is known
for.

###### keywords:

EEG, MEG, source analysis, finite element method, source modeling

[ibb_label] organization=Institute for Biomagnetism and Biosignalanalysis,
University of Münster, city=Münster, country=Germany

[imi_label] organization=Institute of Medical Informatics, University of
Münster, city=Münster, country=Germany

[occ_label] organization=Otto Creutzfeldt Center for Cognitive and Behavioral
Neuroscience, University of Münster, city=Münster, country=Germany

[upjv_label] organization=INSERM U1105, Research Group on Multimodal Analysis
of Brain Function, Jules Verne University of Picardie, city=Amiens,
country=France

[chup_label] organization=Pediatric functional exploration of the nervous
system department, CHU Picardie, city=Amiens, country=France

[tum_label] organization=Institute of History and Ethics in Medicine,
Technical University of Munich, city=Munich, country=Germany

[tuni_label] organization=Computing Sciences Unit, Faculty of Information
Technology and Communication Sciences, Tampere University, city=Tampere,
country=Finland

[umit_label] organization=Institute of Electrical and Biomedical Engineering,
Private University for Health Sciences, Medical Informatics and Technology,
city=Hall in Tyrol, country=Austria

[wwu_label] organization=Faculty of Mathematics and Computer Science,
University of Münster, city=Münster, country=Germany

##  1 Introduction

In electro- (EEG) and magnetoencephalography (MEG), the so-called EEG and MEG
forward problems consist of simulating the sensor measurements that a given
neural activity would generate [1, 2]. The neural activity is typically
modeled as a linear combination of mathematical point dipoles, see e.g. [1, 2,
3]. The problem of forward modeling is thus the computation of electric
potentials (EEG) and magnetic fluxes (MEG) evoked by a single point dipole.
Due to the complex geometry of the head, an accurate approximation of
potentials and fluxes relies on a realistic volume conductor model, where one
then employs some numerical method [1, 2, 4]. There are various approaches to
the computation of these approximations, e.g. boundary element methods (BEM)
[5, 6, 7, 8] or finite difference methods (FDM) [9, 10, 11]. The finite
element method (FEM) has been proposed for modeling complex head geometries
and tissue anisotropy [2, 12, 13] and is therefore in the focus of this work.
A central difficulty for FEM, however, is the embedding of the highly singular
point dipole [14, 15]. We call a strategy for this embedding a FEM potential
approach. Potential approaches can be roughly divided into so-called direct
approaches and subtraction approaches. In direct approaches, the point dipole
is directly incorporated into the FEM right-hand side. State-of-the-art direct
approaches, such as the multipolar Venant approach [16] and the H(div)
approach [17], achieve this by substituting the point dipole with a
regularized dipole-like object. Subtraction approaches on the other hand
analytically handle the dipole by “subtracting” the singularity out of the
problem formulation, arriving at a regular problem, and then post-process the
resulting FEM solution to add the singularity back in. A substantial amount of
work has gone into the development and evaluation of subtraction approaches
[18, 19, 20, 21, 22, 23, 15, 24, 25, 14]. In these studies, subtraction
approaches were shown to posess a rigorous mathematical foundation and to
produce highly accurate results. Despite this, they are typically not used in
practical applications.

The main reason for this is their high computational cost. In source analysis,
one typically needs to compute forward solutions for tens of thousands of
source positions, and, even when employing a fast transfer matrix approach
[26], [15], the time it takes to assemble all FEM right-hand sides for
subtraction-type approaches in realistic head models becomes computationally
infeasible. Different approaches have been suggested to remedy this drawback
of subtraction methods. The projected subtraction approach [23] and the
projected gradient subtraction approach [25] substitute complicated integrals
with simpler ones and thus try to strike a balance between accuracy and
computational burden. The analytic subtraction approach [14] on the other hand
derives analytical formulas for the FEM right-hand side and thus removes the
need for an expensive numerical calculation. While these approaches all reduce
the computational cost of the full subtraction approach as detailed in [15],
none of them has succeeded in making subtraction methods viable in practice.

The fundamental reason for this is that all previous subtraction approaches
lead to FEM formulations that produce dense right-hand sides. Assembling these
right-hand sides requires an iteration over the whole mesh for every dipole
under consideration, where for each element an integration needs to be
performed. Especially for realistic head models, often consisting of millions
of elements (e.g. [27]), this leads to the aforementioned high computational
costs. State-of-the-art direct approaches on the other hand produce sparse
right-hand sides, befitting the focality of a point dipole. These sparse
right-hand sides can then be assembled rapidly.

Aside from the high computational complexity, another problem of the different
subtraction approaches arises for the MEG forward problem. When utilizing a
finite element method to compute magnetic fluxes, one typically first computes
an approximation to the electric potential and then employs the Biot-Savart
law (see [1]). Using a direct approach, this poses no difficulty, as the
numerical solution directly yields the electric current in a form suitable for
the evaluation of Biot-Savart’s law. In contrast, subtraction approaches
produce potentials that are singular at the dipole position, and since Biot-
Savart’s law requires integration over the whole head volume, a direct
application of this law leads to integrals over singular functions, which is
numerically undesirable.

In this paper, we build on top of these observations and propose a new
potential approach, called the localized subtraction approach, that will
alleviate the aforementioned problems. The construction of this approach is
based on the point dipole being only singular in the source position. Hence we
can subtract the singularity out of the problem formulation in a spatially
confined manner, in contrast to previous subtraction approaches, where this
subtraction was performed on the whole head domain. This localization of the
subtraction is the central step in reducing the computational burden of
subtraction approaches.

Furthermore, we present an approach for handling the singular integral arising
in the MEG forward problem when employing a subtraction-type potential
approach. Building on top of the derivation of the classical Geselowitz
formulas [28], the singular part is extracted out of the integral and
rewritten in terms of non-singular surface integrals.

In the methods section, we will derive the corresponding formulas for the EEG
and MEG cases. In particular, we see that for the EEG forward problem the
localized subtraction approach produces right-hand sides that are only nonzero
for degrees of freedom in the close proximity of the dipole, leading to sparse
right-hand sides. This renders the FEM right-hand sides efficiently
computable. In the results section, we will then show that the localized
subtraction approach yields a considerable increase in speed when compared to
previous subtraction approaches. Furthermore, we will show the localized
subtraction approach to be similarly accurate, and in many cases even more
accurate, than other state-of-the-art potential approaches, such as the
analytical subtraction approach [14] and the multipolar Venant approach [16].

##  2 Methods

###  2.1 The EEG forward problem

####  2.1.1 Definition of the problem

Let $\Omega$ be the head domain and $\partial\Omega$ its boundary. In the EEG
forward problem we are interested in determining the electric potential
$u:\Omega\rightarrow\mathbb{R}$, or at least its values at some predetermined
electrode positions, due to some given neural activity. As described in the
introduction, in current practice this reduces to the computation of
potentials generated by mathematical point dipoles. Let
$\sigma:\Omega\rightarrow\mathbb{R}^{3\times 3}$ be the symmetric positive
definite conductivity tensor. We now consider a point dipole at position
$x_{0}\in\Omega$ with moment $M\in\mathbb{R}^{3}$. Using a quasistatic
approximation of Maxwell’s equations [1], the electric potential $u$ due to
the point dipole can be described by

$\displaystyle\operatorname{div}\left(\sigma\nabla u\right)$
$\displaystyle=\operatorname{div}\left(M\cdot\delta_{x_{0}}\right)$ on
$\Omega$, (1) $\displaystyle\langle\sigma\nabla u,\eta\rangle$
$\displaystyle=0$ on $\partial\Omega$, (2)

where $\delta_{x_{0}}$ is the Dirac distribution at $x_{0}$ and $\eta$ is the
unit outer normal of the head domain.

####  2.1.2 Deriving a weak formulation

When trying to solve a partial differential equation using a finite element
approach, the first step consists of deriving a suitable weak formulation. In
doing so, we have to deal with the highly singular term
$\operatorname{div}(M\cdot\delta_{x_{0}})$ arising from (1). The defining
feature of subtraction approaches is that they do not directly incorporate
this term, or an approximation of it, into the weak formulation, but instead
use (1) and (2) to derive a weak formulation for the function $u-\varphi$,
where $\varphi$ is some suitable function which “subtracts” the singularity
out of the problem. We will now discuss how such a $\varphi$ can be
constructed.

We assume the conductivity $\sigma$ to be constant on a neighborhood of
$x_{0}$. Let $\sigma^{\infty}\in\mathbb{R}^{3\times 3}$ be its value on this
neighborhood. Now consider the equation

$\displaystyle\operatorname{div}\left(\sigma^{\infty}\nabla u^{\infty}\right)$
$\displaystyle=\operatorname{div}\left(M\cdot\delta_{x_{0}}\right)$ on
$\mathbb{R}^{3}$, (3)

which is the equivalent of (1) in an unbounded homogeneous conductor. In this
context, it is straightforward to derive analytical formulas for $u^{\infty}$
(e.g. [15, 23]). Since the right-hand side in (3) is the same as in (1), the
idea is to use this $u^{\infty}$ to eliminate the singularity.

One might now simply set $\varphi=u^{\infty}$. All previous subtraction
methods used this approach, and this is what led to their high computational
demand. The reason for this is that $u^{\infty}$ is non-zero almost
everywhere. Hence, weak formulations, and in a second step FEM
discretizations, derived from this approach contain integrals against
functions that are non-zero on a large portion of the head domain [15, 23].
When evaluating these integrals for the FEM basis functions, a large number of
them will be non-zero, leading to dense right-hand sides.

We instead propose to construct $\varphi$ in a different manner. Since
$\operatorname{div}\left(M\cdot\delta_{x_{0}}\right)$ is only singular at the
source position $x_{0}$, we suggest to only use the local distribution of
$u^{\infty}$ around the source position to remove the singularity from the
problem. We make this more precise. Let $\chi:\Omega\rightarrow\mathbb{R}$ be
a function with $\chi=1$ on a neighborhood of the source position $x_{0}$. The
conceptual idea is to choose $\chi$ in such a way that it is only non-zero in
a small region around $x_{0}$. Then $\chi\cdot u^{\infty}$ is also only non-
zero in a small region around $x_{0}$, ultimately leading to a sparse right-
hand side. Furthermore, it shares the local behavior of $u^{\infty}$ close to
the source position and can hence also be used to eliminate the singularity
from the problem formulation. Depending on whether one aims for a (classical)
continuous Galerkin or a discontinuous Galerkin [29] finite element
discretization, one has to impose different regularity constraints on $\chi$.
We focus on the continuous Galerkin case and refer to [30] for the
discontinuous Galerkin case. In the continuous Galerkin case, we further
require $\chi$ to be continuous and to be contained in the Sobolev space
$H^{1}(\Omega)$. We now make the approach $\varphi=\chi\cdot u^{\infty}$, and
define the correction potential $u^{c}$ by

$u^{c}\coloneqq u-\chi\cdot u^{\infty}.$ (4)

We now derive a weak formulation for this correction potential. The central
step is to understand how the multiplication of $u^{\infty}$ with $\chi$
impacts the distributional derivative. This is answered by the following
lemma, which is proven in A.

######  Lemma 1.

Let $\Omega^{\infty}\subset\Omega$ and $\chi\in H^{1}(\Omega)$ be such that
$\chi$ is identical to $1$ on $\Omega^{\infty}$ and $x_{0}\in\Omega^{\infty}$.
Let $\widetilde{\Omega}\subset\Omega\setminus\Omega^{\infty}$ be a region such
that $\partial\Omega^{\infty}$ is a subset of
$\partial\Omega\cup\partial\widetilde{\Omega}$, and $\chi=0$ outside
$\Omega^{\infty}\cup\widetilde{\Omega}$ (see Figure 1). Then for any test
function $\phi\in C^{\infty}_{c}(\Omega)$111By $C^{\infty}_{c}(\Omega)$ we
denote the space of all infinitely differentiable functions with compact
support in $\Omega$. we have in the distributional sense

$\operatorname{div}\left(\sigma^{\infty}\nabla\left(\chi\cdot
u^{\infty}\right)\right)(\phi)=\\\
\operatorname{div}\left(M\cdot\delta_{x_{0}}\right)(\phi)-\int_{\widetilde{\Omega}}\langle\sigma^{\infty}\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS.$ (5)

Here $\eta$ is the unit outer normal of $\Omega^{\infty}$. Note in particular
that we have $\operatorname{div}\left(M\cdot\delta_{x_{0}}\right)$ on the
right-hand side, and hence can use $\chi\cdot u^{\infty}$ to cancel the
singularity in (1).

Figure 1: 2-dim. visualization of a realization of the construction described
in lemma 1. The whole square is $\Omega$. In the yellow region we have
$\chi=1$, in the white region we have $\chi=0$.

We call $\Omega^{\infty}$ the patch. Furthermore, we call $\widetilde{\Omega}$
the transition region, as $\chi$ transitions from 1 to 0 here. We want to
emphasize that our definition of $\Omega^{\infty}$ differs from the one given
in [15, 23]. In particular, we do not assume the conductivity $\sigma$ to be
constant on the whole patch $\Omega^{\infty}$. Instead, we just assume that
$\Omega^{\infty}$ contains a neighborhood of $x_{0}$ where the conductivity is
constant.

Now let $\sigma^{c}\coloneqq\sigma-\sigma^{\infty}$. Multiplying (1) by a test
function, applying partial integration together with the boundary condition
(2), inserting the decomposition

$\sigma\nabla u=\sigma\nabla u^{c}+\sigma^{c}\nabla\left(\chi\cdot
u^{\infty}\right)+\sigma^{\infty}\nabla\left(\chi\cdot u^{\infty}\right)$ (6)

and using lemma 1, we arrive at the following. Define a bilinear form
$a:H^{1}(\Omega)\times H^{1}(\Omega)\rightarrow\mathbb{R}$ by

$a(w,v)=\int_{\Omega}\langle\sigma\nabla w,\nabla v\rangle\,dV$ (7)

and a linear form $l:H^{1}(\Omega)\rightarrow\mathbb{R}$ by

$l(v)=-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla v\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle v\,dS-\int_{\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla v\rangle\,dV,$ (8)

where $\eta$ is the unit outer normal of $\Omega^{\infty}$. Then the
continuous Galerkin localized subtraction approach is given by the problem of
finding $u^{c}\in H^{1}(\Omega)$ such that for all $v\in H^{1}(\Omega)$ we
have $a(u^{c},v)=l(v)$. For a more detailed derivation, we refer to B.

It is easy to see that $l$ vanishes on constant functions. By standard finite
element theory, described e.g. in [23], one can thus see that the set of
solutions of this problem is of the form

$u^{c}_{0}+\mathbb{R}\cdot 1,$ (9)

reflecting the arbitrary choice of a reference electrode. We have thus derived
a weak formulation and have seen that it produces solutions in the expected
form. Also note that the corresponding family of solutions of (1) and (2)
given by

$u^{c}_{0}+\chi\cdot u^{\infty}+\mathbb{R}\cdot 1$ (10)

does not depend on the choice of $\chi,\Omega^{\infty}$ and
$\widetilde{\Omega}$, as is to be expected. Proofs of these statements can be
found in C. Furthermore, we would like to emphasize that the choice $\chi=1$,
$\Omega^{\infty}=\Omega$ and $\widetilde{\Omega}=\emptyset$ leads to the
classical subtraction approach as described in [18, 19, 20, 21, 22, 23, 15,
24, 25, 14]. Hence the classical subtraction approach can be seen as a special
case of the localized subtraction approach.

####  2.1.3 FEM discretization

To use (8) in a computational algorithm, we need a way to construct
$\chi,\Omega^{\infty}$ and $\widetilde{\Omega}$. We want to propose a strategy
for this. Since we want to use a FEM approach, we can assume that a mesh of
the head domain is given. We then construct $\Omega^{\infty}$ and
$\widetilde{\Omega}$ using so called vertex extensions. If we are given a
certain subset $\mathcal{T}$ of mesh elements, we define the vertex extension
$\mathcal{T}^{\prime}$ of $T$ to be the set of all mesh elements that share at
least one vertex with an element from $\mathcal{T}$. Note that
$\mathcal{T}\subseteq\mathcal{T}^{\prime}$. Now let $x_{0}\in\Omega$ be a
given source position. We then take the mesh element containing $x_{0}$ and
perform a certain number of vertex extensions. We take the result of this as
the patch $\Omega^{\infty}$. We then perform an additional vertex extension on
the patch and define $\widetilde{\Omega}$ to be the union over all elements
contained in the vertex extension of the patch that are not contained in the
patch itself. We then choose $\chi$ as a linear combination of FEM basis
functions which takes the value $1$ on every patch vertex and $0$ on the
remaining vertices. Note that in the common cases of piece-wise linear basis
functions on tetrahedra and piece-wise multilinear basis functions on
hexahedra this already uniquely determines $\chi$. The result of this
construction is illustrated in Figure 1. Note that in this figure two vertex
extensions were used to construct the patch $\Omega^{\infty}$.

This construction leaves the question of how many vertex extensions should be
used to construct the patch. From a computational perspective, it is
beneficial to perform as few vertex extensions as possible since larger
patches lead to a more costly integration during the FEM right-hand side
assembly. Furthermore, for larger patches the accuracy of the approximation
becomes more dependent on the underlying mesh, since e.g. at conductivity
jumps $u^{c}=u-\chi\cdot u^{\infty}$ might exhibit an intricate structure,
requiring a high resolution to properly resolve. On the other hand, choosing
small patches might also lead to suboptimal approximations, since in this case
$\chi\cdot u^{\infty}$ is only non-zero in a small vicinity of the source,
which forces the singular behavior of the potential near the dipole onto the
correction potential $u^{c}$. Hence it becomes more difficult to approximate
$u^{c}$ by finite element trial functions, thus degrading the performance of
the approach. In the Results section, we will investigate how we can optimally
balance these effects.

Finally, we want to elaborate on the computation of the FEM right-hand side.
Let $K$ be a mesh element and $F$ a face of a mesh element. When the right-
hand side assembly is performed element-wise, as is typical in the FEM, and
$\varphi$ is a FEM test function, we need to compute integrals of the form

$\displaystyle I_{T}=$
$\displaystyle\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi\rangle\,dV$ (11) $\displaystyle I_{S}=$
$\displaystyle\int_{F}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\varphi\,dS,$ (12) $\displaystyle I_{P}=$
$\displaystyle\int_{K}\langle\sigma^{c}\nabla
u^{\infty},\nabla\varphi\rangle\,dV,$ (13)

We call $I_{T}$ a transition integral, $I_{S}$ a surface integral and $I_{P}$
a patch integral.

In the case of isotropic $\sigma^{\infty}$ and piecewise affine trial
functions on tetrahedral meshes, analytical expressions for the patch and
surface integrals were derived in [14], yielding the analytical subtraction
approach. We slightly adapted the corresponding formulas to arrive at a more
stable implementation. Furthermore, building on top of [31, 32], we extended
the derivation given in [14] to compute an analytical expression for the
transition integral, enabling a completely analytic computation of the
localized subtraction EEG right-hand side. We refer to D for an in-depth
derivation of these formulas.

If the mesh is not tetrahedral, or the conductivity in the neighborhood of the
source is not isotropic, or the trial functions are not affine, we resort to
numerical integration using the quadrature rules implemented in the DUNE
framework [33], on which the DUNEuro implementation is built [34]. A
discussion on how to choose suitable integration orders can be found in E.

###  2.2 The MEG forward problem

####  2.2.1 Definition of the problem

In magnetoencephalography, one measures magnetic fluxes through magnetometers.
To simulate magnetic fluxes, values of the form

$\int_{F}\langle B,\eta\rangle\,dS,$ (14)

need to be approximated, where $B$ is the magnetic field, $F$ is some surface
enclosed by a pickup coil and $\eta$ is the normal to this surface (see e.g.
[1, 2]). In FEM approaches, this is typically done using some quadrature rule
on the integral (14), and then using Biot-Savart’s law to compute the magnetic
field $B$ at the quadrature points [35]. This law states that, for the
electric current $j$ in the volume conductor $\Omega$, the magnetic field at
some point $x$ outside $\Omega$ is given by

$B(x)=\frac{\mu_{0}}{4\pi}\int_{\Omega}j(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y),$
(15)

where $\mu_{0}$ is the magnetic permeability of the vacuum [1, 2]. In
bioelectromagnetism, the current is typically split into the so-called primary
current $j^{P}$ generated by the neural activity, and the passive volume
currents generated by the electric field and hence given by $-\sigma\nabla u$,
where $u$ is the electric potential [1, 2]. Since the magnetic field depends
linearly on the primary current, we can again concentrate on the case of a
point dipole $j^{P}=M\cdot\delta_{x_{0}}$. Inserting the splitting
$j=j^{P}-\sigma\nabla u$ into (15), we get a splitting $B=B^{P}+B^{S}$ into
the so-called primary magnetic field $B^{P}$ and secondary magnetic field
$B^{S}$. Here, the primary field can be efficiently computed using analytical
expressions and hence poses no problem. The difficult part of the MEG forward
problem is thus to compute the integral

$\int_{\Omega}\sigma(y)\nabla u(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y).$ (16)

####  2.2.2 The boundary subtraction approach

To compute the integral in (16), we need the gradient of the electric
potential $u$, which we can e.g. approximate using a FEM approach. With an
increasing distance to the source, the magnitude of the volume current
decreases. Hence the behavior of $\sigma\nabla u$ close to the source position
has a comparatively high contribution in the integral (16), and it is of
interest to ensure that the FEM approximation is accurate in the vicinity of
the dipole. The local behavior of the numerical solution depends on how the
FEM right-hand side was derived from the dipole, which means it depends on the
potential approach. Since the current state-of-the-art direct approaches, such
as the Venant approach, are based on substituting the dipole with a more
easily manageable object, their local behavior around the source position by
construction deviates from that of a true dipole. Subtraction methods in
contrast incorporate the mathematical point dipole into the problem
formulation and can thus be expected to produce highly accurate approximations
of the true current density in the proximity of the dipole. This is
illustrated in Figure 2, where we see that the current computed using the
localized subtraction approach has the expected dipolar pattern even in the
close vicinity of the source position. Hence, when employing a dipolar source
model for the neural activity, subtraction approaches are a natural choice for
computing currents and their resulting magnetic fields.

Figure 2: Visualization of the volume currents as computed by the localized
subtraction approach. The colors of the lines indicate the magnitude of the
current and the colors in the background show the conductivity of the volume
conductor. The numerical simulation was performed in the multilayer sphere
model mesh_init shown in Figure 3 (a), for a radial dipole with a 1 mm
distance to the conductivity jump. Details on this mesh can be found in the
Results section. The current was visualized by generating 5000 points in a
sphere of radius 2.5 mm around the dipole position and applying a Runge-Kutta
method, using the streamtrace filter implemented in the ParaView software
[36].

This again introduces a singularity into the problem, as subtraction
approaches produce approximations of the form $u_{h}=u_{h}^{c}+\chi\cdot
u^{\infty}$, where $u_{h}^{c}$ is the result of the finite element
discretization. Hence we need to compute

$\int_{\Omega}\sigma(y)\nabla\left(u^{c}_{h}+\chi
u^{\infty}\right)(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y).$ (17)

Since $u^{c}_{h}$ is a linear combination of FEM basis funcitons, it poses no
problem for the integration, which can be performed efficiently using a
transfer matrix approach [15, 26, 23]. The function $\nabla\left(\chi\cdot
u^{\infty}\right)$ on the other hand contains a singularity at the source
position, so that numerical integration is undesirable. One way to deal with
this is projecting either $u^{\infty}$ or its gradient into a space of easily
integrable functions, similar to what is done in [23] or [25]. Since
especially in the vicinity of the singularity the function $u^{\infty}$ is
quite different from its projection, this can introduce a significant error,
as was observed in [25] in the EEG context. Under the assumption of isotropic
$\sigma^{\infty}$ we can instead eliminate the singularity from (17)
altogether. For a set $A$, let $1_{A}$ denote its indicator function, which is
1 on $A$ and $0$ everywhere else. Then, if $\Omega^{\infty}$ and
$\widetilde{\Omega}$ are the patch and transition region, by construction, we
have the decomposition

$\sigma\nabla\left(\chi u^{\infty}\right)=\\\ \sigma\nabla\left(\chi
u^{\infty}\right)\cdot 1_{\widetilde{\Omega}}+\sigma^{c}\nabla u^{\infty}\cdot
1_{\Omega^{\infty}}+\sigma^{\infty}\nabla u^{\infty}\cdot
1_{\Omega^{\infty}}.$ (18)

Note that $\sigma\nabla\left(\chi u^{\infty}\right)\cdot
1_{\widetilde{\Omega}}$ and $\sigma^{c}\nabla u^{\infty}\cdot
1_{\Omega^{\infty}}$ do not contain a singularity, since
$1_{\widetilde{\Omega}}$ and $\sigma^{c}$ vanish on an environment of the
singularity. Thus the singularity of the integrand in (17) can be split off
into $\sigma^{\infty}\nabla u^{\infty}\cdot 1_{\Omega^{\infty}}$. But this
expression can be handled by mimicking the derivation of the classical
Geselowitz formula from [28]. This idea was introduced by [37], where it is
used for MEG forward simulations in the context of the full subtraction
approach [15]. More concretely, we have the following.

######  Lemma 2.

Let $\eta$ be the unit outer normal of $\Omega^{\infty}$. We then have

$\int_{\Omega^{\infty}}\sigma^{\infty}\nabla
u^{\infty}(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y)=\\\
\int_{\partial\Omega^{\infty}}\sigma^{\infty}u^{\infty}(y)\cdot\eta(y)\times\frac{x-y}{\|x-y\|^{3}}\,dS(y).$
(19)

A proof of this lemma can be found in [38], where the authors take special
care to properly handle the singularity at the source position.

In total, we thus arrive at the following approach for computing the secondary
magnetic field.

$\displaystyle\begin{split}-\frac{4\pi}{\mu_{0}}B^{S}(x)=&\int_{\Omega}\sigma(y)\nabla
u^{c}(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y)\\\
+&\int_{\Omega^{\infty}}\sigma^{c}(y)\nabla
u^{\infty}(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y)\\\
+&\int_{\partial\Omega^{\infty}}\sigma^{\infty}u^{\infty}(y)\cdot\eta(y)\times\frac{x-y}{\|x-y\|^{3}}\,dS(y)\\\
+&\int_{\widetilde{\Omega}}\sigma(y)\nabla\left(\chi\cdot
u^{\infty}\right)(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y).\end{split}$ (20)

Note in particular that the right-hand side of (20) contains no singularity.
We call this the boundary subtraction approach for the MEG forward problem
since the key step consists in rewriting the singular volume integral as a
boundary integral.

The first summand in (20) can be assembled using a standard MEG transfer
matrix approach. Similar to the EEG case, we call the second summand patch
flux, the third surface flux, and the fourth transition flux. These can be
assembled in an element-wise manner using numerical integration, which raises
the question of how to choose suitable integration orders. A corresponding
discussion can be found in E.

###  2.3 Implementation

All formulas derived in the Methods section were implemented into the DUNEuro
toolbox [34]222https://gitlab.dune-
project.org/duneuro/duneuro/-/tree/feature/localized-assembler, an open-source
C++ toolbox for neuroscience applications built upon the DUNE framework [39,
33].

###  2.4 Experimental setup

For the numerical simulations in this paper, we used a four-layer sphere
model, where the layers represent the brain, cerebrospinal fluid (CSF), skull,
and skin. The model parameters are described in table 1. In multilayer sphere
models one can derive analytical solutions for the EEG [40] and MEG [41]
forward problems, which we compare against our numerical solutions.

Table 1: Multi-layer sphere parameters Compartment |  Conductivity |  Radius
---|---|---
Brain |  0.33 S/m |  0 \- 78 mm
CSF |  1.79 S/m |  78 - 80 mm
Skull |  0.01 S/m |  80 - 86 mm
Skin |  0.43 S/m |  86 - 92 mm

To employ a finite element approach, one needs to construct a mesh. It is well
known that the mesh construction itself can already introduce a bias in favor
of certain potential approaches [42, 43]. More concretely, in [42] it was
found that meshes that concentrate their nodes in the proximity of
conductivity jumps favor subtraction approaches, while meshes with evenly
distributed nodes favor direct approaches, such as the multipolar Venant
approach.

This can be explained by the fact that the subtraction approach does not solve
for the potential itself, but for the difference in the head potential and the
potential in an infinite homogeneous medium. To produce accurate results, the
FEM needs a high resolution in those areas where the function to be
approximated shows a complicated behavior. In particular, in the proximity of
conductivity jumps we expect the head potential to deviate from the potential
in a homogeneous medium. Hence we expect their difference to exhibit a quite
nontrivial behavior at conductivity jumps, which needs a high mesh resolution
to properly capture. In the brain layer, on the other hand, we expect the
homogeneous potential and the head potential to be quite similar, and hence,
as long as the dipole is not too close to a conductivity jump, a comparatively
low brain resolution is required to achieve good numerical accuracy for the
subtraction approach.

The Venant approach, on the other hand, relies on good node positions for the
placement of monopolar loads, and can hence utilize meshes that are more
refined inside the brain compartment.

For the localized subtraction approach, similar considerations as for the
subtraction approach apply. But in contrast to the subtraction approach, which
solves for the whole difference $u-u^{\infty}$, the localized subtraction
approach solves for $u-\chi\cdot u^{\infty}$, which only differs from $u$ in
the proximity of the source. For dipoles with a large distance to the
conductivity jump, all potential approaches generally produce highly accurate
results, and hence the mesh structure is mainly relevant for the sources close
to a conductivity jump. We thus expect the localized subtraction approach to
produce accurate results as long as the mesh has a high resolution at the
conductivity jumps close to the source positions.

These considerations show the difficulty in choosing a ”fair” mesh for
comparing the accuracies of the different potential approaches. We will use
three different meshes, build in such a way as to provide favorable
circumstances for the different potential approaches introduced before (see
Figure 3). First, a tetrahedral mesh consisting of 787k nodes was constructed
using Gmsh [44] (Figure 3 (a)). This mesh has a high concentration of nodes in
the CSF and skull compartments. We then used the UG software [45] to generate
two more meshes from this initial mesh, one by refining the brain compartment
(benefitting the Venant approach, Figure 3 (b)) and one by refining the skin
compartment (benefitting the subtraction approach, Figure 3 (c)). In the
following, we call the mesh in Figure 3 (a) mesh_init, the one in Figure 3 (b)
mesh_brain and the one in Figure 3 (c) mesh_skin.

Figure 3: Circular sections of the meshes used in the numerical tests. (a)
shows the initial mesh with a high concentration of nodes in the CSF and skull
compartment. (b) shows the mesh from (a) after refining the brain compartment.
(c) shows the mesh from (a) after refining the skin compartment. Mesh (a)
consists of ca. 800k nodes, mesh (b) consists of ca. 1.8 million nodes and
mesh (c) consists of ca. 1.3 million nodes.

For the EEG forward problem we used 200 electrode positions on the sphere
surface, and for the MEG forward problem we used 256 coil positions at a
distance of 110 mm to the sphere center. The electrode positions and the coil
positions respectively were chosen to be approximately uniformly distributed
on their corresponding sphere. We always considered dipoles in the brain
compartment, which we grouped by eccentricity, which we computed as the
distance of the dipole from the center of the sphere divided by the radius of
the brain compartment. Based on physiological reasoning, [29] argues that
sources with a distance of 1-2 mm from the CSF compartment are the most
relevant for the generation of neural electromagnetic fields. We will thus
particularly focus on eccentricities inside this range. Furthermore, it is
well known that the numerical simulation of potentials and fluxes becomes more
difficult for sources close to a conductivity jump, as is e.g. noted in [23].
It can now happen, e.g. due to segmentation inaccuracies, that sources are
placed at a distance of below 1 mm to a conductivity jump. To also cover these
extreme cases, we test sources up to an eccentricity of 0.99, corresponding to
a distance of about 0.8 mm to the conductivity jump.

##  3 Results and Discussion

###  3.1 Patch construction

We first want to investigate the number of vertex extensions one should
perform during the construction of the patch for the localized subtraction
approach (see section 2.1.3 and Figure 1). We generated 1000 tangential and
1000 radial dipoles at an eccentricity of $0.99$, each approximately uniformly
distributed on the corresponding sphere. For each of the three meshes shown in
Figure 3, we then performed an increasing number of vertex extensions for each
dipole and computed the corresponding localized subtraction EEG FEM solution.
We also computed the FEM solutions for the analytical subtraction approach
from [14] for comparison. Once the numerical results were computed, the
relative error with respect to the analytical solution at the electrode
positions was calculated. The results for the tangential dipoles and mesh_init
are shown in Figure 4.

Figure 4: Relative error in the EEG case for 1000 tangential dipoles at 0.99
eccentricity for different numbers of vertex extensions during patch
construction, computed using mesh_init (see Figure 3 (a)). The rightmost
yellow boxplot shows the errors when employing the analytical subtraction
approach from [14].

We see that when going from $0$ extensions, meaning the patch consists only of
the element containing the dipole, to two extensions, the error decreases. It
then stays almost constant up to 5 extensions. The error then significantly
increases from 6 to 8 extensions, and then gradually approaches the error of
the analytical subtraction approach. Note that if enough vertex extensions are
performed to cover the whole volume conductor, the error of the localized
subtraction and the analytical subtraction approach would be identical.
Furthermore, radial dipoles show similar behavior.

The initial decrease can be explained by the singular behavior in the
correction potential, which comes from the need to approximate the dipole
potential in the close proximity of the singularity. Furthermore, for
mesh_init, 6-8 extensions are the point where the patches for dipoles at an
eccentricity of $0.99$ start to grow into the skin compartment. As seen in
Figure 3 (a), the skin compartment of this mesh has a comparatively low
resolution. This makes an approximation of $u-\chi\cdot u^{\infty}$, which
might exhibit quite a complicated behavior at conductivity jumps, difficult.
Indeed, when performing the same experiment for mesh_skin, we see that the
error for larger patches significantly decreases, while the error for smaller
patches barely changes (see Figure 2). In contrast, if one performs this
experiment for mesh_brain and compares against Figure 4, the error for larger
patches barely changes, while the error for smaller patches decreases (see
Figure 1). The reason for this is that small patches can exploit the higher
resolution in the brain area, while larger patches are still limited by the
low skin resolution.

Additionally, when comparing the FEM solutions from a small patch with those
from a large patch, we see that the relative difference is high in the skin
elements close to the source and adjacent to the skull. This also explains
why, after the increase in the error at 6-9 extensions, the error continuously
drops, since with growing patches we need to approximate $u-u^{\infty}$, and
not the more complicated $u-\chi\cdot u^{\infty}$, in the skin elements close
to the source and adjacent to the skull.

In total, these experiments illustrate that larger patch sizes make the
approach more dependent on the global mesh structure, while small patch sizes
mainly depend on the local mesh structure in the proximity of the source. In
particular, this is one of the main advantages of the localized subtraction
approach over the analytical subtraction approach.

Importantly, an increase in vertex extensions is always accompanied by an
increase in computational costs. To achieve the best speed and accuracy, we
thus suggest using two vertex extensions for the patch construction.

###  3.2 Comparison to other potential approaches

We evaluate the localized subtraction approach against the analytical
subtraction approach [14]. Furthermore, we want to test the localized
subtraction approach against the multipolar Venant approach [16], as it is a
widely used and well-performing direct approach. In the following, we
abbreviate the localized subtraction approach as ”loc. sub.”, the analytical
subtraction approach as ”ana. sub.” and the multipolar Venant approach as
”mul. Venant”.

####  3.2.1 Accuracy comparison

We first want to compare the accuracies of the different potential approaches.
For this, we computed forward solutions for radial and tangential dipoles at
different, logarithmically scaled eccentricities, ranging from 0.8803 to
0.9900. At each eccentricity, we used 1000 approximately uniformly distributed
dipoles and computed the relative error with respect to the analytic
solutions. For the EEG forward simulation, we shifted all forward solutions to
have zero mean. For the MEG forward simulation, we computed the whole magnetic
field vector. Furthermore, we only solved the MEG forward problem for
tangential dipoles, since radial dipoles produce no magnetic field outside of
a spherically symmetric volume conductor [41]. The results for mesh_init (see
Figure 3 (a)) are shown in Figure 5 and Figure 6.

Figure 5: Accuracy comparison of EEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for radial dipoles at different eccentricities using mesh_init (see Figure 3
(a)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

Figure 6: Accuracy comparison of MEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for tangential dipoles at different eccentricities using mesh_init (see Figure
3 (a)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

We first discuss the EEG case. For radial sources, we see in Figure 5 that in
mesh_init the localized subtraction approach outperforms the analytical
subtraction and multipolar Venant approaches over all investigated
eccentricities and produces a median error of below $1\%$ even for sources at
$0.99$ eccentricity (meaning a distance of $0.8$ mm to the conductivity jump).
Note that we computed the right-hand sides for the analytical subtraction
approach using the analytical formulas from [14], their larger errors are thus
not due to an inaccurate integration but instead result from a relatively low
resolution in the skin compartment, as described in section 3.1. The EEG
results for tangential dipoles show similar behavior, with the errors being
overall smaller (see Figure 1).

Performing the same comparison in mesh_brain shows that, when compared to the
results in Figure 5, refining the brain leads to a large decrease in the
errors of the multipolar Venant approach, a small decrease for the localized
subtraction approach and a barely noticeable decrease for the analytical
subtraction approach (see Figures 2 and 3 for the results and Figure 3 (b) for
the mesh). Furthermore, in this mesh, the localized subtraction approach and
the multipolar Venant approach show approximately the same accuracy. Comparing
Figure 5 against the corresponding results obtained using mesh_skin, we
instead see that the errors for the multipolar Venant approach and the
localized subtraction approach barely change, while the errors for the
analytical subtraction approach notably decrease (see Figures 5 and 6 for the
results and Figure 3 (c) for the mesh). In particular, in this mesh, the
analytical subtraction approach even surpasses the multipolar Venant approach
at some eccentricities with regard to accuracy but does not reach the accuracy
of the localized subtraction approach at any eccentricity.

These results illustrate the limitations of the different potential
approaches. Refining the brain compartment enables a more favorable placement
of monopolar loads, which leads to a decrease in errors for the multipolar
Venant approach when going from mesh_init to mesh_brain. Since the resolution
in the brain compartment is the limiting factor for the placement of monopolar
loads, this also explains why the error for the multipole approach hardly
changes when refining the skin compartment. The analytical subtraction
approach on the other hand shows the opposite behavior. In this case, the low
resolution of the skull-skin interface in mesh_init is the limiting factor,
and hence we see a large error decrease after refining the skin, and only a
very minor decrease after refining the brain. Similar to the analytical
subtraction approach, the localized subtraction approach depends on a high
resolution in the region where the correction potential shows a complicated
behavior. But in contrast to the analytical subtraction approach, this region
is confined to the close proximity of the source positions, and hence the high
CSF resolution in mesh_init already suffices to achieve a good approximation
of $u^{c}=u-\chi\cdot u^{\infty}$. In particular, this explains why the errors
for the localized subtraction approach only slightly decrease when refining
either the brain or the skin.

We now discuss the MEG case. Figure 6 shows the accuracy comparisons for
mesh_init. We see that, in this mesh, both the analytical subtraction and
localized subtraction approach produce highly accurate results and are in
particular more accurate than the multipolar Venant approach in the region of
1-2 mm distance to the conductivity jump. It is especially noteworthy that the
analytical subtraction approach is able to produce such accurate results, as
it produced larger errors in the EEG forward problem in this mesh. The reason
for the EEG errors not having a noticeable influence on the MEG accuracy is
that the errors in the EEG forward problem for the analytical subtraction
approach are mainly located in the skin compartment. But since the skin
compartment is distant from the source, with the low-conducting skull in
between, the volume currents in the skin are low in magnitude and hence their
contribution to the magnetic field is relatively small. In particular, these
results demonstrate the strong influence of the volume currents in the close
proximity of the source on the magnetic field. The accuracy of the analytical
subtraction and localized subtraction MEG forward solutions also shows the
impact of faithful current modeling, which was already suggested by Figure 2.
Furthermore, in contrast to the EEG case, the larger errors of the analytical
subtraction approach for the higher eccentricities are entirely due to the
numerical quadrature. When choosing a sufficiently large integration order for
both approaches, the MEG analytical subtraction approach333 Note that the
analytical subtraction approaches analytically computes the EEG right hand
sides. If one then uses the Biot-Savart law to compute the magnetic field, a
numerical integration is still required. is as accurate as the localized
subtraction approach. In Figure 6, for the analytical subtraction approach an
integration order of $5$ was chosen to achieve a somewhat reasonable
computation time, while the integration orders for the localized subtraction
approach were chosen as described in E.

When performing the same comparison for mesh_skin, we see that the errors are
essentially the same as in Figure 6 (see Figure 7). If we furthermore look at
a single potential approach and compare the errors in mesh_init and
mesh_brain, we see that for all three investigated approaches the errors were
smaller in mesh_brain. Moreover, in mesh_brain the different potential
approaches are comparable with regard to accuracy (see Figure 4).

The results in mesh_skin and mesh_brain again highlight the dominant role of
the volume currents in the close vicinity of the source. Refining the skin
does not enable a better representation of the volume currents close to the
sources while refining the brain does. This explains why the errors do not
change when going from mesh_init to mesh_skin, but decrease when going from
mesh_init to mesh_brain.

Furthermore, these results are in line with previous investigations, where it
was shown that the MEG forward problem is most sensitive to conductivity
uncertainties in the proximity of the dipole location [26], and is insensitive
to the conductivities of the skull and the skin [2].

####  3.2.2 Efficiency comparison

Finally, we now want to compare the computational effort of the different
potential approaches. For this we measured the time it took to solve the EEG
and MEG forward problems for $1000$ right-hand sides for sources at an
eccentricity of $0.99$, assuming the transfer matrix has already been
computed. Hence, we are measuring the time it takes to assemble the right-hand
sides and then compute the product of these right-hand sides with the transfer
matrix. We chose an eccentricity of $0.99$ since the effort needed for the
localized subtraction right-hand side assembly depends on the distance of the
source to the conductivity jump. One needs to assemble the patch, the
boundary, and the transition integrals and fluxes as described in (11) - (13)
and (20). The computational effort for the surface and transition integrals
does not depend on the eccentricity, but the effort of computing the patch
integrals increases for sources closer to the conductivity jump. More
concretely, the patch integrals vanish for every element in the same
compartment as the source, as there we have $\sigma^{c}=0$, and hence we see a
larger number of non-zero patch integrals when increasing the eccentricity, as
a larger portion of the patch grows outside the brain compartment.
Furthermore, for highly eccentric sources the patch elements with a non-zero
integral get close to the source, and hence we need a high integration order
on these elements.444We refer to E for an in-depth discussion of integration
orders. By measuring the computation time for sources with an eccentricity of
$0.99$, which is slightly above the physiologically relevant range, we thus
arrive at a worst-case estimate for the runtime of the localized subtraction
approach. The experiments were run on an AMD Ryzen Threadripper 3960X CPU. The
results for the different meshes (see Figure 3) are shown in the tables 2 \-
4. We see that the time needed for the multipole approach and for the EEG
localized subtraction approach are essentially independent of the mesh
size555Note that the mesh size of course still influences the time needed for
the computation of the transfer matrices.. Furthermore, we see that the
computation time for the MEG localized subtraction approach is the smallest
for mesh_brain. This can be explained by the fact that a higher resolution in
the brain compartment leads to a larger portion of patch elements lying in the
same compartment as the source, which eliminates the need for an expensive
numerical quadrature on these elements. And most importantly, we see that in
all meshes the localized subtraction approach is orders of magnitude faster
than the analytical subtraction approach. Hence the localized subtraction
approach achieves its central goal of being a computationally efficient
approach based on subtraction of the singularity. Furthermore, note that the
time needed for the computation of the transfer matrices is considerably
larger than the time needed for the subsequent localized subtraction right-
hand side assembly. Since for state-of-the-art direct potential approaches,
such as the multipolar Venant approach [16], the total time needed to solve
the forward problems is also dominated by the time needed to compute the
transfer matrices, the localized subtraction approach is thus a viable
alternative to direct FEM potential approaches in practical applications.

Table 2: Forward simulation times for $1000$ sources at an eccentricity of $0.99$ in mesh_init Approach |  loc. sub. |  ana. sub. |  mul. Venant
---|---|---|---
EEG time (s) |  $0.93$ |  $655$ |  $0.02$
MEG time (s) |  $13$ |  $18925$ |  $0.04$
Table 3: Forward simulation times for $1000$ sources at an eccentricity of $0.99$ in mesh_brain Approach |  loc. sub. |  ana. sub. |  mul. Venant
---|---|---|---
EEG time (s) |  $0.96$ |  $1396$ |  $0.02$
MEG time (s) |  $6.8$ |  $43175$ |  $0.04$
Table 4: Forward simulation times for $1000$ sources at an eccentricity of $0.99$ in mesh_skin Approach |  loc. sub. |  ana. sub. |  mul. Venant
---|---|---|---
EEG time (s) |  $0.98$ |  $1142$ |  $0.02$
MEG time (s) |  $12.8$ |  $30709$ |  $0.04$

##  4 Conclusion

We presented a new potential approach based on a localization of the well-
known subtraction approach. We investigated the theory of this new approach
and validated and compared it against established approaches. In spherical
head models, the new approach was shown to rival, and in some cases even
surpass, existing approaches in terms of accuracy while being largely more
efficient than the subtraction approach. We believe the localized subtraction
to be the superior way of analytically dealing with the singularity of the
mathematical point dipole, and a competitive alternative in practical
applications.

## Data/code availability statement

The code implementing the presented approach is available at
https://gitlab.dune-project.org/duneuro/duneuro/-/tree/feature/localized-
assembler. The code used to test different integration orders is available at
https://github.com/MalteHoel/local_integration_test. The data used in the
experiments is available at https://doi.org/10.5281/zenodo.7642826.

## Declaration of competing interests

Declarations of interest: none.

## Credit authorship contribution statement

Malte Höltershinken: Methodology, Software, Formal analysis, Investigation,
Writing - original draft, Writing - review and editing, Visualization. Pia
Lange: Methodology, Software, Investigation, Writing - review and editing. Tim
Erdbrügger: Writing - review and editing. Yvonne Buschermöhle: Writing -
review and editing. Fabrice Wallois: Writing - review and editing, Funding
Acquisition. Alena Buyx: Funding Acquisition. Sampsa Pursiainen: Writing -
review and editing, Funding Acquisition. Johannes Vorwerk: Investigation,
Writing - review and editing, Funding Acquisition. Christian Engwer:
Conceptualization, Methodology, Software, Writing - review and editing,
Supervision. Carsten Wolters: Conceptualization, Methodology, Writing - review
and editing, Supervision, Funding Acquisition.

## Acknowledgements

This work was supported by the Bundesministerium für Gesundheit (BMG) as
project ZMI1-2521FSB006, under the frame of ERA PerMed as project
ERAPERMED2020-227 (MBH, PL, TE, YB, FW, AB, CHW), by the Deutsche
Forschungsgemeinschaft (DFG), project WO1425/10-1 (MBH, TE, YB, CHW), by the
Deutsche Forschungsgemeinschaft through the Cluster of Excellence “Mathematics
Münster: Dynamics - Geometry - Structure” (EXC 2044-390685587) (CE), by the
Academy of Finland (decision 344712) under the frame of ERA PerMed as project
ERAPERMED2020-227 (SP), by the AoF Center of Excellence in Inverse Modelling
and Imaging 2018-2025 (decision 353089) (SP), by the Austrian
Wissenschaftsfonds (FWF), project P 35949 (JV), and by DAAD project 57663920
(MBH, TE, SP, CHW). We acknowledge support from the Open Access Publication
Fund of the University of Muenster.

##  Appendix A Proof of lemma 1

In the following, we assume $\Omega\subset\mathbb{R}^{3}$ to be open.
Furthermore, we denote by $C^{\infty}_{c}(\Omega)$ the space of all infinitely
differentiable, compactly supported functions on $\Omega$. Note that compactly
supported functions on $\Omega$ necessarily vanish near the boundary of
$\Omega$, which implies that boundary integrals involving such functions
vanish. This in particular applies to the boundary integrals arising during
partial integration.

To give a mathematically rigorous derivation of the method, we will make use
of distributions. A distribution (or generalized function) is defined as a
linear functional $T:C^{\infty}_{c}(\Omega)\rightarrow\mathbb{R}$ which is
continuous with respect to a certain topology. Every locally integrable
function $f$ defines a distribution $T_{f}$ via

$T_{f}(\phi)=\int_{\Omega}f\phi\,dV.$

In the following, we will simply write $T_{f}(\phi)=f(\phi)$. A central
property of distributions is that they possess arbitrary derivatives, while
ordinary functions may fail to be (weakly) differentiable. It is thus natural
to state and investigate partial differential equations in a distributional
framework. Nevertheless, we want to emphasize that the following derivation of
the localized subtraction approach does not assume the reader to be familiar
with this theory. For readers interested in an in-depth introduction to
distributions and their application to differential equations, we refer to
[46], chapters 6-8.

For reference, we will first recall lemma 1.

######  Lemma.

Let $\Omega^{\infty}\subset\Omega$ and $\chi\in H^{1}(\Omega)$ be such that
$\chi$ is identical to $1$ on $\Omega^{\infty}$ and $x_{0}\in\Omega^{\infty}$.
Let $\widetilde{\Omega}\subset\Omega\setminus\Omega^{\infty}$ be a region such
that $\partial\Omega^{\infty}$ is a subset of
$\partial\Omega\cup\partial\widetilde{\Omega}$, and $\chi=0$ outside
$\Omega^{\infty}\cup\widetilde{\Omega}$. Then for any test function $\phi\in
C^{\infty}_{c}(\Omega)$ we have

$\operatorname{div}\left(\sigma^{\infty}\nabla\left(\chi\cdot
u^{\infty}\right)\right)(\phi)=\\\
\operatorname{div}\left(M\cdot\delta_{x_{0}}\right)(\phi)-\int_{\widetilde{\Omega}}\langle\sigma^{\infty}\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS.$ (21)

Here $\eta$ is the unit outer normal of $\Omega^{\infty}$.

###### Proof.

First, note that $\chi\cdot u^{\infty}$ is locally integrable, and hence
defines a distribution. We then have by definition

$\displaystyle\operatorname{div}(\sigma^{\infty}\nabla(\chi\cdot
u^{\infty}))(\phi)\stackrel{{\scriptstyle}}{{=}}\int_{\Omega}\chi\cdot
u^{\infty}\operatorname{div}(\sigma^{\infty}\nabla\phi)\,dV$
$\displaystyle=\int_{\Omega}u^{\infty}\operatorname{div}(\sigma^{\infty}\nabla\phi)\,dV+\int_{\Omega}(\chi-1)u^{\infty}\operatorname{div}(\sigma^{\infty}\nabla\phi)\,dV$
$\displaystyle=\operatorname{div}(\sigma^{\infty}\nabla
u^{\infty})(\phi)+\int_{\Omega\setminus\Omega^{\infty}}(\chi-1)u^{\infty}\operatorname{div}(\sigma^{\infty}\nabla\phi)\,dV$
$\displaystyle=\operatorname{div}(M\cdot\delta_{x_{0}})(\phi)+\int_{\Omega\setminus\Omega^{\infty}}(\chi-1)u^{\infty}\operatorname{div}(\sigma^{\infty}\nabla\phi)\,dV.$

Here we have used for the first summand that $u^{\infty}$ solves
$\operatorname{div}(\sigma^{\infty}\nabla
u^{\infty})=\operatorname{div}(M\cdot\delta_{x_{0}})$ on $\mathbb{R}^{3}$ and
that $\phi$ has compact support in $\Omega$. For the second summand, we have
used $\chi=1$ on $\Omega^{\infty}$. Note that $(\chi-1)u^{\infty}$ is weakly
differentiable on $\Omega\setminus\Omega^{\infty}$, and hence we can apply
partial integration to the second summand. Letting $\eta$ denote the
corresponding unit outer normal, this yields

$\int_{\Omega\setminus\Omega^{\infty}}(\chi-1)u^{\infty}\operatorname{div}(\sigma^{\infty}\nabla\phi)\,dV\\\
=\int_{\partial(\Omega\setminus\Omega^{\infty})}(\chi-1)u^{\infty}\langle\sigma^{\infty}\nabla\phi,\eta\rangle\,dS\\\
-\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{\infty}\nabla((\chi-1)u^{\infty}),\nabla\phi\rangle\,dV$

Since
$\partial(\Omega\setminus\Omega^{\infty})\subset\partial\Omega\cup\partial\Omega^{\infty}$
and we have $\nabla\phi=0$ on $\partial\Omega$ and $\chi-1=0$ on
$\partial\Omega^{\infty}$ the boundary integral vanishes. For the volume
integral, we get by again applying partial integration

$-\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{\infty}\nabla((\chi-1)u^{\infty}),\nabla\phi\rangle\,dV\\\
=-\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{\infty}\nabla(\chi\cdot
u^{\infty}),\nabla\phi\rangle\,dV+\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\nabla\phi\rangle\,dV\\\
=-\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{\infty}\nabla(\chi\cdot
u^{\infty}),\nabla\phi\rangle\,dV+\int_{\partial(\Omega\setminus\Omega^{\infty})}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS\\\
-\int_{\Omega\setminus\Omega^{\infty}}\operatorname{div}(\sigma^{\infty}\nabla
u^{\infty})\phi\,dV.$

Since on $\mathbb{R}^{3}\setminus\\{x_{0}\\}$ we have
$\operatorname{div}(\sigma^{\infty}\nabla u^{\infty})=0$, the last integral
vanishes. For the boundary integral, since $\phi$ has compact support in
$\Omega$ the portion of the integral on the boundary of $\Omega$ vanishes.
Noting that the outer normal of $\Omega\setminus\Omega^{\infty}$ on the
boundary of $\Omega^{\infty}$ in the interior of $\Omega$ is the negative of
the outer normal of $\Omega^{\infty}$, we arrive at (21). We want to emphasize
that this derivation also covers the case where the patch touches the boundary
of $\Omega$ since in this case, the corresponding portion of the boundary
integral in (21) vanishes since $\phi$ has compact support. This completes the
proof. ∎

##  Appendix B Derivation of the weak formulation

Let

$a(w,v)=\int_{\Omega}\langle\sigma\nabla w,\nabla v\rangle\,dV$ (22)

and

$l(v)=-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla v\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle v\,dS-\int_{\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla v\rangle\,dV$ (23)

be defined as in the main paper. Then $a$ defines a continuous bilinear form
and $l$ defines a continuous linear form on $H^{1}(\Omega)$. We want to show
that, for $u^{c}\in H^{1}(\Omega)$ a solution of $a(u^{c},v)=l(v)$ for all
$v\in H^{1}(\Omega)$, the function $u^{c}+\chi\cdot u^{\infty}$ solves the
original equation in a suitable sense. Reading this backward can serve as a
derivation of the weak formulation.

######  Lemma B1.

Let $u^{c}\in H^{1}(\Omega)$ be a solution of $a(u^{c},v)=l(v)$ for all $v\in
H^{1}(\Omega)$. Then $u=u^{c}+\chi\cdot u^{\infty}$ solves

$\displaystyle\operatorname{div}(\sigma\nabla u)$
$\displaystyle=\operatorname{div}(M\cdot\delta_{x_{0}})$ on $\Omega$
$\displaystyle\langle\sigma\nabla u,\eta\rangle$ $\displaystyle=0$
$\displaystyle\text{ on $\partial\Omega$}.$

###### Proof.

Since $u^{c}\in H^{1}(\Omega)$, $u^{\infty}$ is locally integrable and $\chi$
is bounded, it follows that $u$ is locally integrable and hence defines a
distribution. We then have for a test function $\phi\in
C^{\infty}_{c}(\Omega)$

$\operatorname{div}(\sigma\nabla u)(\phi)=\operatorname{div}(\sigma\nabla
u^{c})(\phi)+\operatorname{div}(\sigma^{c}\nabla(\chi\cdot
u^{\infty}))(\phi)\\\ +\operatorname{div}(\sigma^{\infty}\nabla(\chi\cdot
u^{\infty}))(\phi)\\\ =-\int_{\Omega}\langle\sigma\nabla
u^{c},\nabla\phi\rangle\,dV-\int_{\Omega}\langle\sigma^{c}\nabla(\chi\cdot
u^{\infty}),\nabla\phi\rangle\,dV\\\
+\operatorname{div}(M\cdot\delta_{x_{0}})(\phi)-\int_{\widetilde{\Omega}}\langle\sigma^{\infty}\nabla(\chi\cdot
u^{\infty}),\nabla\phi\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS\\\
=\operatorname{div}(M\cdot\delta_{x_{0}})(\phi)-a(u^{c},\phi)-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS-\int_{\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla\phi\rangle\,dV\\\
=\operatorname{div}(M\cdot\delta_{x_{0}})(\phi),$

where the last equality follows from $a(u^{c},\phi)=l(\phi)$ and the
definition of $l$. We have thus shown that $u$ solves

$\displaystyle\operatorname{div}(\sigma\nabla u)$
$\displaystyle=\operatorname{div}(M\cdot\delta_{x_{0}})$ $\displaystyle\text{
on $\Omega$}.$

Now note that since $\operatorname{div}(\sigma^{\infty}\nabla u^{\infty})=0$
on $\mathbb{R}^{3}\setminus\\{x_{0}\\}$, we have for a smooth function
$\phi\in C^{\infty}\left(\overline{\Omega}\right)$ with
$\operatorname{supp}(\phi)\cap\\{x_{0}\\}=\emptyset$ that

$\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS=\int_{\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\nabla\phi\rangle\,dV.$

We thus have

$l(\phi)=-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\phi\,dS-\int_{\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla\phi\rangle\,dV\\\
=-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV-\int_{\Omega^{\infty}}\langle\sigma\nabla
u^{\infty},\nabla\phi\rangle\,dV\\\
=-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV-\int_{\Omega^{\infty}}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV\\\
=-\int_{\Omega}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\phi\rangle\,dV.$

Inserting this into $a(u^{c},\phi)-l(\phi)=0$ yields

$\int_{\Omega}\langle\sigma\nabla u,\nabla\phi\rangle\,dV=0.$

Since $\operatorname{div}(\sigma\nabla u)=0$ on $\Omega\setminus\\{x_{0}\\}$
applying partial integration again yields

$\int_{\partial\Omega}\langle\sigma\nabla u,\eta\rangle\phi\,dS=0,$

which yields $\langle\sigma\nabla u,\eta\rangle=0$ on $\partial\Omega$. This
completes the proof. ∎

##  Appendix C Well-posedness of the problem

In the main paper, we claimed that the functional $l$ vanishes on constant
functions and that the family of solutions to the total potential $u$ does not
depend on the choice of $\chi,\Omega^{\infty}$, and $\widetilde{\Omega}$. We
want to give proof of these claims.

######  Lemma C1.

The functional $l$ in (23) vanishes on constant functions.

###### Proof.

Using $\nabla 1=0$ we get

$l(1)=0\iff\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS=0.$

Now let $R>0$ be large enough so that $\Omega^{\infty}$ is contained inside
the interior of the sphere $S_{R}(x_{0})$ of radius $R$ around $x_{0}$. Let
$W\subset\mathbb{R}^{3}$ denote region bounded by $\Omega^{\infty}$ and
$S_{R}(x_{0})$. Since $\operatorname{div}(\sigma^{\infty}\nabla u^{\infty})=0$
on $\mathbb{R}^{3}\setminus\\{x_{0}\\}$ the divergence theorem then implies

$0=\int_{W}\operatorname{div}\left(\sigma^{\infty}\nabla
u^{\infty}\right)\,dV=\int_{\partial W}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS,$

where $\eta$ denotes the outer normal of $W$. Since we have $\partial W$ =
$S_{R}(x_{0})\cup\partial\Omega^{\infty}$, and the outer normal of $W$ on
$\partial\Omega^{\infty}$ is the negative of the outer normal of
$\Omega^{\infty}$, this implies

$\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS=\int_{S_{R}(x_{0})}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS,$

where $\eta$ again denotes the corresponding unit outer normals. We now have
for $x\in\mathbb{R}^{3}\setminus\\{x_{0}\\}$ that

$\displaystyle
u^{\infty}(x)=\frac{1}{4\pi\sqrt{\det\sigma^{\infty}}}\cdot\frac{\langle
M,\left(\sigma^{\infty}\right)^{-1}(x-x_{0})\rangle}{\|\left(\sigma^{\infty}\right)^{-\frac{1}{2}}(x-x_{0})\|^{3}},$

see e.g. [23]. Denoting $R=\|x-x_{0}\|$, a straightforward computation shows
that $\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\in\mathcal{O}\left(\frac{1}{R^{3}}\right)$. This
implies

$\int_{S_{R}(x_{0})}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS\in\mathcal{O}\left(\frac{1}{R}\right),$

and hence

$\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS=\int_{S_{R}(x_{0})}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS\stackrel{{\scriptstyle
R\rightarrow\infty}}{{\longrightarrow}}0.$

This is only possible if

$\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS=0.$

∎

Note that the previous lemma essentially consists of showing that

$\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\,dS=0.$

The same integral also appears in [23], where the authors use a symmetry
argument to deduce that it is equal to $0$. Since our argument in lemma C1
takes a somewhat different approach, we deemed it interesting enough to
include.

######  Lemma C2.

Let $\Omega^{\infty}_{1},\widetilde{\Omega}_{1},\chi_{1}$ and
$\Omega^{\infty}_{2},\widetilde{\Omega}_{2},\chi_{2}$ be two choices of patch,
transition region, and cutoff function. Let $u^{c}_{1}$ and $u^{c}_{2}$ be
solutions to the corresponding weak formulations. Then we have

$u^{c}_{1}+\chi_{1}\cdot u^{\infty}+\mathbb{R}\cdot 1=u^{c}_{2}+\chi_{2}\cdot
u^{\infty}+\mathbb{R}\cdot 1.$

The potentials defined by the localized subtraction approach thus do not
depend on the choice of the patch, transition region, and cutoff function.

###### Proof.

It obviously suffices to show that any localized subtraction approach produces
the same solution family as the ordinary subtraction approach. We can thus
without loss of generality assume $\Omega^{\infty}_{2}=\Omega$,
$\widetilde{\Omega}_{2}=\emptyset$ and $\chi_{2}=1$.

Now note that

$u^{c}_{1}-u^{c}_{2}+(\chi_{1}-1)\cdot u^{\infty}\in H^{1}(\Omega),$

since $\chi_{1}-1$ vanishes on an environment of $x_{0}$. Then we have for
$v\in H^{1}(\Omega)$

$a(u^{c}_{1}-u^{c}_{2}+(\chi_{1}-1)\cdot u^{\infty},v)\\\
=a(u^{c}_{1},v)-a(u^{c}_{2},v)+\int_{\Omega}\langle\sigma\nabla(\chi_{1}-1)u^{\infty},\nabla
v\rangle\,dV.$ (24)

Since $u^{c}_{1}$ and $u^{c}_{2}$ are solutions to their respective weak
problems, we get

$a(u^{c}_{1},v)-a(u^{c}_{2},v)=-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi_{1}\cdot
u^{\infty}\right),\nabla v\rangle\,dV\\\
-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle v\,dS-\int_{\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla v\rangle\,dV\\\ +\int_{\Omega}\langle\sigma^{c}\nabla
u^{\infty},\nabla
v\rangle\,dV+\int_{\partial\Omega}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle v\,dS\\\
=\left(\int_{\partial\Omega}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle
v\,dS-\int_{\partial\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle v\,dS\right)\\\
+\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla
v\rangle\,dV-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi_{1}\cdot
u^{\infty}\right),\nabla v\rangle\,dV\\\
=\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{\infty}\nabla
u^{\infty},\nabla
v\rangle\,dV+\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma^{c}\nabla
u^{\infty},\nabla v\rangle\,dV\\\
-\int_{\widetilde{\Omega}}\langle\sigma\nabla\left(\chi_{1}\cdot
u^{\infty}\right),\nabla v\rangle\,dV\\\
=\int_{\Omega\setminus\Omega^{\infty}}\langle\sigma\nabla\left(1-\chi_{1})u^{\infty}\right),\nabla
v\rangle\,dV.$

Inserting this into (24) and noting that $\chi_{1}-1=0$ on $\Omega^{\infty}$,
we see that

$a(u_{1}^{c}-u_{2}^{c}+(\chi_{1}-1)\cdot u^{\infty},v)=0.$

This implies that $(u_{1}^{c}+\chi_{1}\cdot
u^{\infty})-(u_{2}^{c}+u^{\infty})=c$ for some constant $c\in\mathbb{R}$,
completing the proof. ∎

##  Appendix D Analytical expressions for EEG FEM

In the main paper, we noted that in the case of affine trial functions on
tetrahedral meshes we compute the FEM right-hand side using analytical
expression. We have to assemble three types of integrals, namely

$\displaystyle I_{P}=$ $\displaystyle\int_{K}\langle\sigma^{c}\nabla
u^{\infty},\nabla\varphi\rangle\,dV,$ (25) $\displaystyle I_{S}=$
$\displaystyle\int_{F}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\varphi\,dS,$ (26) $\displaystyle I_{T}=$
$\displaystyle\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi\rangle\,dV,$ (27)

where $K$ is some tetrahedron and $F$ is some triangle. We call $I_{P}$ a
patch integral, $I_{S}$ a surface integral and $I_{T}$ a transition integral.
One can find analytical expressions for the patch integral and the surface
integral in [14]. We slightly modified those expressions and used the
techniques described in [14, 31, 32] to derive an analytical expression for
the transition integral777If the reader is mainly interested in implementing
the formulas for an analytic computation of the transition integral, all
necessary information is given in the following. But since this work is an
extension of [14, 31, 32], we strongly recommend first studying these papers..
To set the stage for these formulas we first have to introduce some notation.
In the following, we always assume a dipolar source at position $x_{0}$ with
moment $M$, and we assume $x_{0}\not\in F$.

Let $F$ be some triangle. Let $p_{1},p_{2},p_{3}$ be an arbitrary numbering of
the vertices of $F$. Then set

$\displaystyle s_{1}=$
$\displaystyle\frac{p_{3}-p_{2}}{\|p_{3}-p_{2}\|},\gamma_{1}=\|p_{3}-p_{2}\|,$
$\displaystyle s_{2}=$
$\displaystyle\frac{p_{1}-p_{3}}{\|p_{1}-p_{3}\|},\gamma_{2}=\|p_{1}-p_{3}\|,$
$\displaystyle s_{3}=$
$\displaystyle\frac{p_{2}-p_{1}}{\|p_{2}-p_{1}\|},\gamma_{3}=\|p_{2}-p_{1}\|.$

Now set $u=s_{3},w=\frac{s_{1}\times s_{2}}{\|s_{1}\times s_{2}\|}$ and
$v=w\times u$, and let

$m_{i}=s_{i}\times w$

for $i\in\\{1,2,3\\}$. Then set

$\displaystyle u_{3}=$ $\displaystyle\langle u,p_{3}-p_{1}\rangle,$
$\displaystyle v_{3}=$ $\displaystyle\langle v,p_{3}-p_{1}\rangle,$
$\displaystyle u_{0}=$ $\displaystyle\langle u,x_{0}-p_{1}\rangle,$
$\displaystyle v_{0}=$ $\displaystyle\langle v,x_{0}-p_{1}\rangle,$
$\displaystyle w_{0}=$ $\displaystyle-\langle w,x_{0}-p_{1}\rangle.$

Then define

$\displaystyle t_{1}=$
$\displaystyle\frac{v_{0}(u_{3}-\gamma_{3})+v_{3}(\gamma_{3}-u_{0})}{\gamma_{1}},$
$\displaystyle t_{2}=$
$\displaystyle\frac{u_{0}v_{3}-v_{0}u_{3}}{\gamma_{2}},$ $\displaystyle
t_{3}=$ $\displaystyle v_{0}.$

and

$\displaystyle\gamma_{1}^{-}=$
$\displaystyle-\frac{(\gamma_{3}-u_{0})(\gamma_{3}-u_{3})+v_{0}v_{3}}{\gamma_{1}}$
$\displaystyle\gamma_{1}^{+}=$
$\displaystyle\frac{(u_{3}-u_{0})(u_{3}-\gamma_{3})+v_{3}(v_{3}-v_{0})}{\gamma_{1}}$
$\displaystyle\gamma_{2}^{-}=$
$\displaystyle-\frac{u_{3}(u_{3}-u_{0})+v_{3}(v_{3}-v_{0})}{\gamma_{2}}$
$\displaystyle\gamma_{2}^{+}=$
$\displaystyle\frac{u_{0}u_{3}+v_{0}v_{3}}{\gamma_{2}}$
$\displaystyle\gamma_{3}^{-}=$ $\displaystyle-u_{0}$
$\displaystyle\gamma_{3}^{+}=$ $\displaystyle\gamma_{3}-u_{0}.$

We refer to [14] for a geometric intuition behind these values (note that
there is an error in the formula for $\gamma_{2}^{-}$ in [14], compare e.g.
against [32]). Essentially, we have constructed a map of the triangle and
computed the coordinates of its vertices and parametrizations of its edges. We
then set for $1\leq i\leq 3$

$\displaystyle R_{i}=$ $\displaystyle\sqrt{t_{i}^{2}+w_{0}^{2}}$
$\displaystyle R_{i}^{\pm}=$
$\displaystyle\sqrt{t_{i}^{2}+w_{0}^{2}+\left(\gamma_{i}^{\pm}\right)^{2}}$
$\displaystyle f_{i}=$
$\displaystyle\begin{cases}\ln\left(\frac{R_{i}^{+}+\gamma_{i}^{+}}{R_{i}^{-}+\gamma_{i}^{-}}\right)&\text{if
$\gamma_{i}^{-}\geq 0$}\\\
\ln\left(\frac{R_{i}^{-}-\gamma_{i}^{-}}{R_{i}^{+}-\gamma_{i}^{+}}\right)&\text{if
$\gamma_{i}^{-}<0$}\end{cases}$ (28) $\displaystyle\beta_{i}=$
$\displaystyle\arctan\left(\frac{t_{i}\gamma_{i}^{+}}{R_{i}^{2}+|w_{0}|R_{i}^{+}}\right)$
$\displaystyle-\arctan\left(\frac{t_{i}\gamma_{i}^{-}}{R_{i}^{2}+|w_{0}|R_{i}^{-}}\right)$
$\displaystyle\beta=$ $\displaystyle\beta_{1}+\beta_{2}+\beta_{3}.$

Note in particular that our definition of $f_{i}$ (see equation (28)) differs
from the one given in [14, 31, 32], where no case distinction is made and the
first case is always taken. Carefully looking at the formulas, we see that
$R_{i}^{-}$ is the distance of the source position $x_{0}$ to the point
$p_{i+1}$ and $R_{i}^{+}$ is the distance from the point $p_{i+2}$ (with
$p_{4}:=p_{1},p_{5}:=p_{2}$). Furthermore, $\gamma_{i}^{-}$ is the oriented
distance from the projection of the source position onto the affine line
defined by $p_{i+1}$ and $p_{i+2}$, where the positive orientation is given by
$p_{i+2}-p_{i+1}$, to the point $p_{i+1}$. Now assume that the source position
lies (approximately) on a line defined by some edge of the triangle and
happens to lie on the ray defined by the edge and the positive direction. In
this case, we have $R_{i}^{+}\approx-\gamma_{i}^{+}$ and
$R_{i}^{-}\approx-\gamma_{i}^{-}$ for this edge, and the computation of
$f_{i}$ using the first case becomes unstable or even undefined. Taking a
close look at the derivation of the formulas in [14, 31, 32], we see that
$f_{i}$ arises from the computation of the integral

$\int_{\gamma_{i}^{-}}^{\gamma_{i}^{+}}\frac{1}{\sqrt{w_{0}^{2}+t_{i}^{2}+t^{2}}}\,dt.$

A straightforward computation reveals that for $a>0$ we have

$\frac{d}{dx}\left(\frac{1}{2}\ln\left(\frac{\sqrt{x^{2}+a}+x}{\sqrt{x^{2}+a}-x}\right)\right)=\frac{1}{\sqrt{x^{2}+a}},$

which can be used for the computation of the aforementioned integral.
Expanding the argument of the logarithm by
$\frac{\sqrt{x^{2}+a}+x}{\sqrt{x^{2}+a}+x}$ leads to the first expression for
$f_{i}$ given in equation (28), while expanding with
$\frac{\sqrt{x^{2}+a}-x}{\sqrt{x^{2}+a}-x}$ leads to the second expression. If
$w_{0}\neq 0$ or $t_{i}\neq 0$, the computation of $f_{i}$ is unproblematic.
Otherwise the assumption $x_{0}\not\in F$ forces either
$\gamma_{i}^{+},\gamma_{i}^{-}>0$ or $\gamma_{i}^{+},\gamma_{i}^{-}<0$. In the
first case, the first expression for $f_{i}$ provides a stable way to compute
$f_{i}$, and in the second case, the second expression for $f_{i}$ provides a
stable way to compute $f_{i}$. The formula for $f_{i}$ given above thus avoids
problems arising from the case where the source position lies on the line
defined by some edge of the triangle.

Now set

$a=\begin{pmatrix}-\frac{1}{\gamma_{3}}\\\ \frac{1}{\gamma_{3}}\\\
0\end{pmatrix},b=\begin{pmatrix}\frac{\frac{u_{3}}{\gamma_{3}}-1}{v_{3}}\\\
-\frac{u_{3}}{\gamma_{3}v_{3}}\\\
\frac{1}{v_{3}}\end{pmatrix},\varphi(u_{0},v_{0})=\begin{pmatrix}1\\\ 0\\\
0\end{pmatrix}+u_{0}a+v_{0}b.$

We can now give the analytical formulas for the transition integrals (see
equation (27)).

######  Lemma D1.

We want to compute the integral

$\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi\rangle\,dV,$

where $K$ is some tetrahedron, $\varphi$ is a nodal Lagrange basis function,
and $\chi$ is affine on $K$. Let $F_{1},F_{2},F_{3},F_{4}$ be the faces of the
tetrahedron. Let $1\leq j\leq 4$. Let $p_{1}(F_{j}),p_{2}(F_{J}),p_{3}(F_{j})$
be some numbering of the vertices of $F_{j}$. Then set

$c(F_{j})=\begin{pmatrix}\chi(p_{1}(F_{j}))\\\ \chi(p_{2}(F_{j}))\\\
\chi(p_{3}(F_{j}))\end{pmatrix}.$

Furthermore, let $M(F_{j})=(M_{0},M_{u},M_{v})\in\mathbb{R}^{3\times 3}$,
where

$\displaystyle M_{0}=$ $\displaystyle\operatorname{sign}(w_{0})\beta
w-\sum_{i=1}^{3}f_{i}m_{i},$ $\displaystyle M_{u}=$
$\displaystyle\sum_{i=1}^{3}\langle
u,s_{i}\rangle\left(f_{i}t_{i}s_{i}-(R_{i}^{+}-R_{i}^{-})m_{i}\right)$
$\displaystyle-|w_{0}|\beta u-\left(\sum_{i=1}^{3}\langle u,m_{i}\rangle
f_{i}\right)w_{0}w,$ $\displaystyle M_{v}=$
$\displaystyle\sum_{i=1}^{3}\langle
v,s_{i}\rangle\left(f_{i}t_{i}s_{i}-(R_{i}^{+}-R_{i}^{-})m_{i}\right)$
$\displaystyle-|w_{0}|\beta v-\left(\sum_{i=1}^{3}\langle u,m_{i}\rangle
f_{i}\right)w_{0}w,$

where it is understood that all values in the above equations depend on the
face $F_{j}$ and are computed as defined above. Then set

$I(F_{j})=M(F_{j})\cdot\begin{pmatrix}\langle\varphi(u_{0},v_{0}),c(F_{j})\rangle\\\
\langle a,c(F_{j})\rangle\\\ \langle b,c(F_{j})\rangle\end{pmatrix}.$

Now let $\mu_{K}:\widehat{K}\rightarrow K$ be the affine map from the
reference tetrahedron to $K$. Let $A\in\mathbb{R}^{3\times 3}$ be its
(constant) jacobian. Let
$\widehat{\varphi}_{1},\widehat{\varphi}_{2},\widehat{\varphi}_{3},\widehat{\varphi}_{4}$
be the shape functions on the reference element and let
$G=(\nabla\widehat{\varphi}_{1},\nabla\widehat{\varphi}_{2},\nabla\widehat{\varphi}_{3},\nabla\widehat{\varphi}_{4})\in\mathbb{R}^{3\times
4}$ be the matrix containing their (constant) jacobians. Denote by
$\eta(F_{j})$ the (constant) unit outer normal of the face $F_{j}$.

Then we have

$\begin{pmatrix}\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{1}\rangle\,dV\\\
\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{2}\rangle\,dV\\\
\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{3}\rangle\,dV\\\
\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{4}\rangle\,dV\end{pmatrix}\\\
=\frac{\left(\sigma
A^{-T}G\right)^{T}}{4\pi\sigma^{\infty}}\sum_{j=1}^{4}\langle
I(F_{j}),M\rangle\eta(F_{j}).$

###### Proof.

We have

$\begin{pmatrix}\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{1}\rangle\,dV\\\
\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{2}\rangle\,dV\\\
\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{3}\rangle\,dV\\\
\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi_{4}\rangle\,dV\end{pmatrix}\\\
=\begin{pmatrix}\nabla\varphi_{1}^{T}\\\ \nabla\varphi_{2}^{T}\\\
\nabla\varphi_{3}^{T}\\\
\nabla\varphi_{4}^{T}\end{pmatrix}\cdot\sigma\cdot\int_{K}\nabla\left(\chi\cdot
u^{\infty}\right)\,dV\\\ =\frac{\left(\sigma
A^{-T}G\right)^{T}}{4\pi\sigma^{\infty}}\int_{K}\nabla\left(\chi\cdot\frac{\langle
M,x-x_{0}\rangle}{\|x-x_{0}\|^{3}}\right)\,dV,$

where in the last line we have used the chain rule, the symmetry of $\sigma$,
and the definition of $u^{\infty}$. Now the divergence theorem yields

$\int_{K}\nabla\left(\chi\cdot\frac{\langle
M,x-x_{0}\rangle}{\|x-x_{0}\|^{3}}\right)\,dV\\\
=\sum_{j=1}^{4}\int_{F_{j}}\chi\cdot\frac{\langle
M,x-x_{0}\rangle}{\|x-x_{0}\|^{3}}\cdot\eta(F_{j})\,dS\\\
=\sum_{j=1}^{4}\eta(F_{j})\langle
M,\int_{F_{j}}\chi\cdot\frac{x-x_{0}}{\|x-x_{0}\|^{3}}\,dS\rangle.$

We thus need to compute integrals of the form

$\int_{F}\chi\cdot\frac{x-x_{0}}{\|x-x_{0}\|^{3}}\,dS.$

for a triangle $F$. Let $p_{1},p_{2},p_{3}$ denote the vertices of this
triangle and let $\varphi_{i}$ be the unique affine function that is $1$ on
$p_{i}$ and zero on the remaining vertices. Since by construction $\chi$ is
affine on the face $F$, we have

$\chi=\sum_{i=1}^{3}\chi(p_{i})\cdot\varphi_{i}=\langle\begin{pmatrix}\varphi_{1}\\\
\varphi_{2}\\\ \varphi_{3}\end{pmatrix},c(F)\rangle.$

Now set $O=(u,v,w)$ and define

$\Phi:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3};(z_{1},z_{2})\mapsto
p_{1}+O\cdot\begin{pmatrix}z_{1}\\\ z_{2}\\\ 0\end{pmatrix}.$

Let
$\widehat{F}^{\prime}:=\operatorname{conv}\\{(0,0),(\gamma_{3},0),(u_{3},v_{3})\\}$.
The values defined above were chosen in such a way that
$\Phi|_{\widehat{F}^{\prime}}$ defines a map of the triangle $F$. To make the
following computation easier we follow [14] and introduce an additional
translation

$\tau_{\rho}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3};z\mapsto z+\rho,$

where $\rho=(u_{0},v_{0})$ are the coordinates of the projection of the source
position $x_{0}$ onto the plane defined by the triangle. We now set
$\widehat{F}=\tau_{\rho}^{-1}(\widehat{F}^{\prime})$ and set
$\Phi_{\rho}=\Phi\circ\tau_{\rho}$. Then $\Phi_{\rho}|_{\widehat{F}}$ defines
a map of $F$. Note that the gramian of this map is $1$, and hence we can
compute

$\int_{F}\chi\cdot\frac{x-x_{0}}{\|x-x_{0}\|^{3}}\,dS.=\int_{\widehat{F}}\chi(\Phi_{\rho}(z))\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz.$

One can easily check that

$\begin{pmatrix}\varphi_{1}\\\ \varphi_{2}\\\
\varphi_{3}\end{pmatrix}(\Phi(z_{1},z_{2}))=(e_{1},a,b)\cdot\begin{pmatrix}1\\\
z_{1}\\\ z_{2}\end{pmatrix},$

where $e_{1}=(1,0,0)$ is the first unit vector, see e.g. [14], and hence

$\begin{pmatrix}\varphi_{1}\\\ \varphi_{2}\\\
\varphi_{3}\end{pmatrix}(\Phi_{\rho}(z_{1},z_{2}))=\varphi(u_{0},v_{0})+z_{1}a+z_{2}b,$

where $\varphi(u_{0},v_{0})$ is the value defined just before the statement of
the lemma. Putting this together, we get

$\int_{F}\chi\cdot\frac{x-x_{0}}{\|x-x_{0}\|^{3}}\,dS=\int_{\widehat{F}}\chi(\Phi_{\rho}(z))\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz\\\
=\int_{\widehat{F}}\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz\cdot\langle\varphi(u_{0},v_{0}),c(F)\rangle\\\
+\int_{\widehat{F}}z_{1}\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz\cdot\langle
a,c(F)\rangle\\\
+\int_{\widehat{F}}z_{2}\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz\cdot\langle
b,c(F)\rangle.$

Hence we need to compute the values

$\displaystyle M_{0}$
$\displaystyle=\int_{\widehat{F}}\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz$
$\displaystyle M_{u}$
$\displaystyle=\int_{\widehat{F}}z_{1}\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz$
$\displaystyle M_{v}$
$\displaystyle=\int_{\widehat{F}}z_{2}\frac{\Phi_{\rho}(z)-x_{0}}{\|\Phi_{\rho}(z)-x_{0}\|^{3}}\,dz.$

The value of $M_{0}$ reported in the statement of the lemma is then a direct
consequence of combining (25), (26), and (30) from [32] (note that we followed
the convention of [14], meaning our definition of $w_{0}$ and the definition
of [32] differ by a factor of $-1$). Similarly, the values of $M_{u}$ and
$M_{v}$ can be computed by combining (25), (27), (31), and (32) from [32].
This completes the derivation. ∎

##  Appendix E Choosing integration orders

In the main paper it was described that if we want to solve the EEG forward
problem and are given a tetrahedral mesh, piecewise affine test functions, and
isotropic conductivity in the element containing the source position, we can
assemble the localized subtraction FEM right-hand side analytically. If these
conditions for the EEG forward problem are not met, e.g. if we use a
hexahedral mesh, or if we want to solve the MEG forward problem, we resort to
numerical integration. This of course raises the question of choosing suitable
integration orders, which we want to discuss in this section.

Let $K$ be a mesh element and let $F$ be a face of a mesh element. In the EEG
case, we are concerned with the integrals

$\displaystyle I_{P}=$ $\displaystyle\int_{K}\langle\sigma^{c}\nabla
u^{\infty},\nabla\varphi\rangle\,dV,$ (29) $\displaystyle I_{S}=$
$\displaystyle\int_{F}\langle\sigma^{\infty}\nabla
u^{\infty},\eta\rangle\varphi\,dS,$ (30) $\displaystyle I_{T}=$
$\displaystyle\int_{K}\langle\sigma\nabla\left(\chi\cdot
u^{\infty}\right),\nabla\varphi\rangle\,dV,$ (31)

where $\varphi$ is a test function. We call $I_{P}$ a patch integral, $I_{S}$
a surface integral and $I_{T}$ a transition integral. In the MEG case, we are
concerned with the integrals

$\displaystyle F_{P}=$ $\displaystyle\int_{K}\sigma^{c}(y)\nabla
u^{\infty}(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y),$ (32) $\displaystyle
F_{S}=$
$\displaystyle\int_{F}\sigma^{\infty}u^{\infty}(y)\cdot\eta(y)\times\frac{x-y}{\|x-y\|^{3}}\,dS(y),$
(33) $\displaystyle F_{T}=$
$\displaystyle\int_{K}\sigma(y)\nabla\left(\chi\cdot
u^{\infty}\right)(y)\times\frac{x-y}{\|x-y\|^{3}}\,dV(y),,$ (34)

where we call $F_{P}$ a patch flux, $F_{S}$ a surface flux and $F_{T}$ a
transition flux.

In [14], the author investigated the error introduced by numerical quadrature
for the patch integrals and surface integrals. In particular, the author shows
experimentally that the central quantity for assessing the accuracy of the
numerical quadrature is the ratio $d/a$, where $d$ is the distance of the
source position to the element under consideration and $a$ is the edge length
of this element. Here, smaller ratios demand higher integration orders to
achieve the same accuracy.

This observation is indeed a mathematical necessity. When performing numerical
quadrature, one typically chooses a map of the element $K$ (resp. the face
$F$) to pull back the integral to the reference element, and then performs the
numerical integration on the reference element. A straightforward computation
now shows that when scaling the source position and the element $K$ (resp. the
face $F$) with a scalar $\lambda\neq 0$, the integrands (29)-(31), when pulled
back to the reference element, are scaled by a factor of $1/\lambda$, at least
for affine and multilinear test functions. This implies that the analytic and
numeric integrals are also scaled by a factor of $1/\lambda$. This in turn
implies that the relative error introduced by numerical integration is scale-
invariant. If the coil position is also scaled, another straightforward
computation shows that the integrands (32)-(34) on the reference element are
scaled by a factor of $1/\lambda^{2}$, again implying scale-invariance of the
relative error introduced by numerical integration.

This observation enables us to specify rules for choosing integration orders
independently of a concrete mesh since the scale-invariance of the relative
error enables us to investigate the accuracy of different integration orders
at a reference scale.

Our tests showed that, when compared to the dipole distance or the maximum
edge length, the element shape and the dipole orientation only had a minor
influence on the relative integration error. Among the elements of fixed
maximum edge length, the equilateral ones tended to produce the largest
errors, which we suppose is a consequence of them having the lowest density of
quadrature points888The underlying heuristic here is that a higher density of
quadrature points should enable a more accurate numerical integration. The
density of quadrature points is given by the number of quadrature points
divided by the element volume. When performing numerical integration by
pulling back to a reference element, the number of quadrature points for a
given integration order is independent of the element shape, and hence the
quadrature point density for a given integration order only depends on the
reciprocal of the element volume.. We will thus derive quadrature rules by
investigating the relative error introduced by numerical integration for a
single equilateral element of side length $a=1$ for different dipoles at
distance $d=d/a$.

More concretely, we used equilateral tetrahedral and hexahedral elements and
computed the relative error introduced by numerical integration for different
integration orders and different values of $\frac{d}{a}$. In particular, we
used the tetrahedron with vertices

$\displaystyle p_{1}$
$\displaystyle=\begin{pmatrix}0,&-\frac{1}{2},&-\frac{\sqrt{3}}{6}\end{pmatrix},$
$\displaystyle p_{2}$
$\displaystyle=\begin{pmatrix}0,&\frac{1}{2},&-\frac{\sqrt{3}}{6}\end{pmatrix},$
$\displaystyle p_{3}$
$\displaystyle=\begin{pmatrix}0,&0,&-\frac{\sqrt{3}}{3}\end{pmatrix},$
$\displaystyle p_{4}$
$\displaystyle=\begin{pmatrix}-\frac{\sqrt{6}}{3},&0,&0\end{pmatrix},$

and for the surface integrals the triangle with vertices $p_{1},p_{2},p_{3}$.

Furthermore, we took the cube with vertices

$\displaystyle q_{1}$
$\displaystyle=\begin{pmatrix}-1,&-0.5,-0.5\end{pmatrix}$ $\displaystyle
q_{2}$ $\displaystyle=\begin{pmatrix}0,&-0.5,&-0.5\end{pmatrix}$
$\displaystyle q_{3}$
$\displaystyle=\begin{pmatrix}-1,&0.5,&-0.5\end{pmatrix}$ $\displaystyle
q_{4}$ $\displaystyle=\begin{pmatrix}0,&0.5,&-0.5\end{pmatrix}$ $\displaystyle
q_{5}$ $\displaystyle=\begin{pmatrix}-1,&-0.5,&0.5\end{pmatrix}$
$\displaystyle q_{6}$ $\displaystyle=\begin{pmatrix}0,&-0.5,&0.5\end{pmatrix}$
$\displaystyle q_{7}$ $\displaystyle=\begin{pmatrix}-1,&0.5,&0.5\end{pmatrix}$
$\displaystyle q_{8}$ $\displaystyle=\begin{pmatrix}0,&0.5,&0.5\end{pmatrix}$

and for the surface integrals the square possessing the vertices
$q_{2},q_{4},q_{6},q_{8}$. We used the source position

$x_{0}=d\cdot\begin{pmatrix}1,&0.1,&0.1\end{pmatrix}.$

Note that we did not choose $(1,0,0)$, since for this position some
integration orders perform unusually well due to what we suppose are symmetry
reasons. Furthermore, we used the point

$\begin{pmatrix}20,0,0\end{pmatrix}$

as a coil position for the computation of the magnetic field. Our experiments
suggest that, as long as the distance of the source to the element is
significantly smaller than the distance of the coil to the element, the
concrete coil position only has a minor influence on the error. Furthermore,
with increasing distance of the coil to the element, the error tended to
decrease. Considering the setup of modern MEG machines, it thus seems
justified to use a single coil at a distance of $2$cm from the element in our
experiments. Relative errors were computed by comparing against pseudo-
analytic integrals obtained by using an integration order of 50. For the
transition integrals and fluxes, we defined $\chi$ to be the affine resp. the
multilinear function which is 1 on all vertices belonging to the face closest
to the dipole and 0 on the remaining vertices. The code used for the
computations in this section was implemented as a DUNE
module999https://github.com/MalteHoel/local_integration_test.

Since we perform two vertex extensions, we can expect the distance-edge length
ratios for the transition and surface integrals to be distributed around $2$.
To check our intuition, we took mesh_init from the results section of the main
paper and $1000$ dipoles at $0.99$ eccentricity. For each of these dipoles, we
built the corresponding patch and computed for every element in the transition
region the ratio d/a. The result is shown in Figure 1. We computed $a$ as the
longest edge length and $d$ as the minimum value of the distances of the
source position to the corners and the centers of the faces of the
tetrahedron.

Figure 1: Histogram visualizing the relative frequency of distance-edge length
ratios during the computation of surface and transition integrals ((33) and
(34)) in a tetrahedral mesh. To generate this image, $1000$ dipoles at an
eccentricity of $0.99$ were considered.

We see that the distribution peaks around $2$, and only a very small fraction
of element-source pairs have a ratio less than $1$, as we expected. Hence if
we choose an integration order that is accurate for a ratio of $1$ (and
slightly below $1$), we can expect the relative error for the surface and
transition integrals to be small.

We first want to derive integration orders for the MEG fluxes (32)-(34) in the
case of a tetrahedral mesh. In Figures 2 and 3 we see the relative errors for
different integration orders for the surface (33) and transition fluxes (34).
Based on this, we suggest using an integration order of $6$ for the surface
integrals and an integration order of $5$ for the transition integrals.

Figure 2: Relative error introduced by numerical integration for the MEG
surface flux (33) for a tetrahedral element for different integration orders.
On the $x$-axis we have the distance-edge length ratio and on the $y$-axis we
have the relative error.

Figure 3: Relative error introduced by numerical integration of the MEG
transition flux (34) for a tetrahedral element for different integration
orders. On the $x$-axis we have the distance-edge length ratio and on the
$y$-axis we have the relative error

For the patch fluxes, we cannot expect d/a to be nicely bounded below.
Instead, we suggest choosing an individual integration order for each element-
source combination. More concretely, once d/a has been computed, one can use
the results shown in Figure 4 to choose an appropriate integration order. This
leads to the rule already described in table 1. This rule gives accurate
integration orders as long as $\frac{d}{a}\geq\frac{1}{6}$, which should
suffice for all meshes and sources used in practical applications.

Figure 4: Relative error introduced by numerical integration of the MEG patch flux for (32) a tetrahedral element for different integration orders. On the $x$-axis we have the distance-edge length ratio and on the $y$-axis we have the relative error Table 1: Tetrahedral MEG patch flux integration order $\rule{0.0pt}{8.5359pt}\frac{\text{distance}}{\text{edge length}}$ $\left(\frac{d}{a}\right)$ |  $\frac{d}{a}\geq 0.5$ |  $\frac{d}{a}\geq 0.4$ |  $\frac{d}{a}\geq 0.33$ |  $\frac{d}{a}\geq 0.25$ |  $\frac{d}{a}\geq 0.17$
---|---|---|---|---|---
Intorder |  8 |  9 |  11 |  13 |  20

We now consider the hexahedral case. In this case, we need to specify
integration orders for the EEG and the MEG case. We again perform two vertex
extensions to construct the patch. Note that we can expect a larger lower
bound for the $d/a$-ratio for the surface and transition integrals since
hexahedral meshes are typically generated from voxel-based segmentation
followed by a potential geometry adaption. To check this assumption, we
constructed a geometry-adapted 1 mm hexahedral mesh with a node shift factor
of $0.33$, as described in [47]. Similar to the tetrahedral case, we then took
1000 dipoles at an eccentricity of $0.99$ and constructed the corresponding
patches and transition regions. For each dipole and each element inside the
corresponding transition region, we then computed the ratio $d/a$, using the
same approach as in the tetrahedral case. The result is shown in Figure 5. We
see that the distribution peaks at $2.5$, and that the vast majority of
elements have a ratio above $1.5$. Moreover, not a single element had a ratio
below $1$. If we thus choose integration orders which are accurate up to a
ratio of $1$, we can expect the numerical integration to only produce a small
relative error for the surface and transition integrals (30), (31) and fluxes
(33), (34).

Furthermore, it was discussed in the main paper that physiologically relevant
sources have a distance of at least $1$ mm to the conductivity jump. Under the
assumption that hexahedral meshes used in practice have an edge length of at
most $2$ mm, an integration order that is accurate up to a ratio of $0.5$ thus
suffices for an accurate computation of the patch integral (29) and patch flux
(32).101010Note however that we strongly recommend using a 1 mm mesh over a 2
mm mesh.

Figure 5: Histogram visualizing the relative frequency of distance-edge length
ratios during the computation of surface and transition integrals ((33) and
(34)) in a geometry-adapted hexahedral mesh. To generate this image, $1000$
dipoles at an eccentricity of $0.99$ were considered.

Similar to the tetrahedral case, we computed the relative error introduced by
numerical integration on a single hexahedral element for different
$d/a$-values. Figures 6, 7, and 8 show the corresponding results for the EEG
case, while Figures 9, 10, and 11 show the results for the MEG case.

Based on the discussion in the previous paragraphs and the results in these
figures, we thus suggest choosing an integration order according to table 2.
Note that we can be conservative when choosing orders since when using the
localized subtraction approach the computation times are small enough to not
be a serious concern anymore.

Figure 6: Relative error introduced by numerical integration of the EEG
surface integral (30) for a hexahedral element for different integration
orders. On the $x$-axis we have the distance-edge length ratio and on the
$y$-axis we have the relative error

Figure 7: Relative error introduced by numerical integration of the EEG
transition integral (31) for a hexahedral element for different integration
orders. On the $x$-axis we have the distance-edge length ratio and on the
$y$-axis we have the relative error

Figure 8: Relative error introduced by numerical integration of the EEG patch
integral (29) for a hexahedral element for different integration orders. On
the $x$-axis we have the distance-edge length ratio and on the $y$-axis we
have the relative error

Figure 9: Relative error introduced by numerical integration of the MEG
surface flux (33) for a hexahedral element for different integration orders.
On the $x$-axis we have the distance-edge length ratio and on the $y$-axis we
have the relative error

Figure 10: Relative error introduced by numerical integration of the MEG
transition flux (34) for a hexahedral element for different integration
orders. On the $x$-axis we have the distance-edge length ratio and on the
$y$-axis we have the relative error

Figure 11: Relative error introduced by numerical integration of the MEG patch flux (32) for a hexahedral element for different integration orders. On the $x$-axis we have the distance-edge length ratio and on the $y$-axis we have the relative error Table 2: Hexahedral integration orders |  Surface |  Transition |  Patch
---|---|---|---
EEG |  6 |  6 |  8
MEG |  6 |  6 |  8

##  Appendix F Vertex Extension Comparison

In the main paper, we investigated the influence of the number of vertex
extensions on the accuracy of the localized subtraction approach. In
particular, we referred to this appendix for the figures showing the
corresponding error curves for mesh_brain and mesh_skin. The results for
mesh_brain can be seen in Figure 1, and the results for mesh_skin can be seen
in Figure 2.

Figure 1: Relative error in the EEG case for 1000 tangential dipoles at 0.99
eccentricity for different numbers of vertex extensions during patch
construction, computed using mesh_brain (see Figure 3 (b)). The rightmost
yellow boxplot shows the errors when employing the analytical subtraction
approach from [14].

Figure 2: Relative error in the EEG case for 1000 tangential dipoles at 0.99
eccentricity for different numbers of vertex extensions during patch
construction, computed using mesh_skin (see Figure 3 (c)). The rightmost
yellow boxplot shows the errors when employing the analytical subtraction
approach from [14].

We also performed the vertex extension comparison for the MEG forward problem.
The corresponding results for mesh_init can be seen in Figure 3. We see an
increase in accuracy when going from 0 to 2 extensions. After that, the errors
barely change. Similar to the EEG case, the initial error decrease is due to
smaller patches forcing more singular behavior onto the correction potential.
But in contrast to the EEG case, the MEG errors do not increase again as the
patch becomes larger. This again demonstrates that the secondary magnetic
field due to a dipolar source is predominantly influenced by the volume
currents in the close proximity of the source. Performing the corresponding
experiment for mesh_brain and mesh_skin leads to the same qualitative results.

Figure 3: Relative error in the MEG case for 1000 tangential dipoles at 0.99
eccentricity for different numbers of vertex extensions during patch
construction, computed using mesh_init (see Figure 3 (a)).

##  Appendix G Accuracy comparison for different potential approaches

In the main paper, we showed the accuracy comparison of the localized
subtraction approach against the multipolar Venant and analytical subtraction
approaches in the EEG case for radial sources and in the MEG cases for
tangential sources in mesh_init. To keep the presentation focused, we moved
the figures showing the remaining comparisons into this appendix. The
comparison for tangential sources in the EEG case in mesh_init can be found in
Figure 1. Figures 2, 3 and 4 show the results for mesh_brain. Figures 5, 6 and
7 show the results for mesh_skin.

Figure 1: Accuracy comparison of EEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for tangential dipoles at different eccentricities using mesh_init (see Figure
3 (a)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

Figure 2: Accuracy comparison of EEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for radial dipoles at different eccentricities using mesh_brain (see Figure 3
(b)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

Figure 3: Accuracy comparison of EEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for tangential dipoles at different eccentricities using mesh_brain (see
Figure 3 (b)). The $y$-axis shows the relative error. The physiologically
relevant sources at 1-2 mm distance from the CSF are highlighted.

Figure 4: Accuracy comparison of MEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for tangential dipoles at different eccentricities using mesh_brain (see
Figure 3 (b)). The $y$-axis shows the relative error. The physiologically
relevant sources at 1-2 mm distance from the CSF are highlighted.

Figure 5: Accuracy comparison of EEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for radial dipoles at different eccentricities using mesh_skin (see Figure 3
(c)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

Figure 6: Accuracy comparison of EEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for tangential dipoles at different eccentricities using mesh_skin (see Figure
3 (c)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

Figure 7: Accuracy comparison of MEG forward simulations using the analytical
subtraction, multipolar Venant, and localized subtraction potential approaches
for tangential dipoles at different eccentricities using mesh_skin (see Figure
3 (c)). The $y$-axis shows the relative error. The physiologically relevant
sources at 1-2 mm distance from the CSF are highlighted.

## References

  * [1] M. Hämäläinen, R. Hari, R. J. Ilmoniemi, J. Knuutila, and O. V. Lounasmaa, “Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain,” Rev. Mod. Phys., vol. 65, pp. 413–497, Apr 1993.
  * [2] R. Brette and A. Destexhe, “Handbook of neural activity measurement,” Cambridge University Press, 2012.
  * [3] J. de Munck, B. van Dijk, and H. Spekreijse, “Mathematical dipoles are adequate to describe realistic generators of human brain activity,” IEEE Transactions on Biomedical Engineering, vol. 35, no. 11, pp. 960–966, 1988\.
  * [4] C. Ramon, P. Schimpf, and J. Haueisen, “Influence of head models on EEG simulations and inverse source localizations,” Biomedical engineering online, vol. 5, p. 10, 02 2006.
  * [5] M. Fuchs, R. Drenckhahn, H. Wischmann, and M. Wagner, “An improved boundary element method for realistic volume-conductor modeling,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 8, pp. 980–997, 1998.
  * [6] S. N. Makarov, M. Hämäläinen, Y. Okada, G. M. Noetscher, J. Ahveninen, and A. Nummenmaa, “Boundary element fast multipole method for enhanced modeling of neurophysiological recordings,” IEEE Transactions on Biomedical Engineering, vol. 68, no. 1, pp. 308–318, 2021.
  * [7] 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.
  * [8] Y. Ataseven, Z. Akalin Acar, C. Acar, and N. Gencer, “Parallel implementation of the accelerated BEM approach for EMSI of the human brain,” Medical and biological engineering and computing, vol. 46, pp. 671–9, 08 2008\.
  * [9] 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, pp. 95–111, 01 2014.
  * [10] S. Turovets, V. Volkov, A. Zherdetsky, A. Prakonina, and A. D. Malony, “A 3D finite-difference BiCG iterative solver with the fourier-jacobi preconditioner for the anisotropic EIT/EEG forward problem,” Computational and Mathematical Methods in Medicine, vol. 2014, 2014.
  * [11] E. Cuartas, C. Acosta-Medina, G. Castellanos-Dominguez, and D. Mantini, “A finite-difference solution for the EEG forward problem in inhomogeneous anisotropic media,” Brain Topography, vol. 32, 03 2019.
  * [12] H. Azizollahi, M. Darbas, M. M. Diallo, A. E. Badia, and S. Lohrengel, “EEG in neonates: Forward modeling and sensitivity analysis with respect to variations of the conductivity,” Mathematical Biosciences and Engineering, vol. 15, no. 4, pp. 905–932, 2018.
  * [13] M. Vermaas, M. Piastra, T. Oostendorp, N. Ramsey, and P. Tiesinga, “FEMfuns: A volume conduction modeling pipeline that includes resistive, capacitive or dispersive tissue and electrodes,” Neuroinformatics, vol. 18, 10 2020.
  * [14] L. Beltrachini, “The analytical subtraction approach for solving the forward problem in EEG,” Journal of Neural Engineering, vol. 16, p. 056029, sep 2019.
  * [15] F. Drechsler, C. Wolters, T. Dierkes, H. Si, and L. Grasedyck, “A full subtraction approach for finite element method based source analysis using constrained delaunay tetrahedralisation,” NeuroImage, vol. 46, pp. 1055–65, 04 2009.
  * [16] J. Vorwerk, A. Hanrath, C. H. Wolters, and L. Grasedyck, “The multipole approach for EEG forward modeling using the finite element method,” NeuroImage, vol. 201, p. 116039, 2019.
  * [17] T. Miinalainen, A. Rezaei, D. us, A. Nüßing, C. Engwer, C. Wolters, and S. Pursiainen, “A realistic, accurate and fast source modeling approach for the EEG forward problem,” NeuroImage, vol. 184, 08 2018.
  * [18] M. Thevenet, O. Bertrand, F. Perrin, T. Dumont, and J. Pernier, “The finite element method for a realistic head model of electrical brain activities: preliminary results,” Clinical Physics and Physiological Measurement, vol. 12, p. 89, jan 1991.
  * [19] M. P. SP van den Broek, H Zhou, “Computation of neuromagnetic fields using finite-element method and Biot-Savart law,” Meg. Biol. Eng. Comput., vol. 34, pp. 21–26, jan 1996.
  * [20] K. Awada, D. Jackson, J. Williams, D. Wilton, S. Baumann, and A. Papanicolaou, “Computational aspects of finite element modeling in EEG source localization,” IEEE Transactions on Biomedical Engineering, vol. 44, no. 8, pp. 736–752, 1997.
  * [21] G. Marin, C. Guerin, S. Baillet, L. Garnero, and G. Meunier, “Influence of skull anisotropy for the forward and inverse problem in EEG: Simulation studies using FEM on realistic head models,” Human Brain Mapping, vol. 6, no. 4, pp. 250–269, 1998.
  * [22] P. Schimpf, C. Ramon, and J. Haueisen, “Dipole models for the EEG and MEG,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 5, pp. 409–418, 2002.
  * [23] C. Wolters, H. Köstler, C. Möller, J. Härdtlein, L. Grasedyck, and W. Hackbusch, “Numerical mathematics of the subtraction method for the modeling of a current dipole in EEG source reconstruction using finite element head models,” SIAM J. Scientific Computing, vol. 30, pp. 24–45, 01 2007.
  * [24] C. Engwer, J. Vorwerk, J. Ludewig, and C. H. Wolters, “A discontinuous Galerkin method to solve the EEG forward problem using the subtraction approach,” SIAM Journal on Scientific Computing, vol. 39, no. 1, pp. B138–B164, 2017.
  * [25] L. Beltrachini, “Sensitivity of the projected subtraction approach to mesh degeneracies and its impact on the forward problem in EEG,” IEEE Transactions on Biomedical Engineering, vol. 66, no. 1, pp. 273–282, 2019.
  * [26] N. Gencer and C. Acar, “Sensitivity of EEG and MEG measurements to tissue conductivity,” Physics in medicine and biology, vol. 49, pp. 701–17, 04 2004.
  * [27] M. C. Piastra, S. Schrader, A. Nüßing, M. Antonakakis, T. Medani, A. Wollbrink, C. Engwer, and C. H. Wolters, “The WWU DUNEuro reference data set for combined EEG/MEG source analysis,” June 2020.
  * [28] D. Geselowitz, “On the magnetic field generated outside an inhomogeneous volume conductor by internal current sources,” IEEE Transactions on Magnetics, vol. 6, no. 2, pp. 346–347, 1970.
  * [29] M. C. Piastra, A. Nüßing, J. Vorwerk, H. Bornfleth, R. Oostenveld, C. Engwer, and C. H. Wolters, “The discontinuous Galerkin finite element method for solving the MEG and the combined MEG/EEG forward problem,” Frontiers in Neuroscience, vol. 12, 2018.
  * [30] A. Nuessing, Fitted and unfitted finite element methods for solving the EEG forward problem.  dissertation, University of Muenster, 2018.
  * [31] D. Wilton, S. Rao, A. Glisson, D. Schaubert, O. Al-Bundak, and C. Butler, “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Transactions on Antennas and Propagation, vol. 32, no. 3, pp. 276–281, 1984.
  * [32] R. 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\.
  * [33] P. Bastian, M. Blatt, A. Dedner, N.-A. Dreier, C. Engwer, R. Fritze, C. Gräser, C. Grüninger, D. Kempf, R. Klöfkorn, M. Ohlberger, and O. Sander, “The Dune framework: Basic concepts and recent developments,” Computers & Mathematics with Applications, vol. 81, pp. 75–112, 2021.  Development and Application of Open-source Software for Problems with Numerical PDEs.
  * [34] S. Schrader, A. Westhoff, M. C. Piastra, T. Miinalainen, S. Pursiainen, J. Vorwerk, H. Brinck, C. H. Wolters, and C. Engwer, “Duneuro—a software toolbox for forward modeling in bioelectromagnetism,” PLOS ONE, vol. 16, pp. 1–21, 06 2021.
  * [35] B. Roth and S. Sato, “Accurate and efficient formulas for averaging the magnetic field over a circular coil,” Biomagnetism: Clinical aspects, jan 1992.
  * [36] J. P. Ahrens, B. Geveci, and C. C. Law, “ParaView: An end-user tool for large-data visualization,” in The Visualization Handbook, 2005.
  * [37] M. C. Piastra, New Finite Element Methods for Solving the MEG and the Combined MEG/EEG Forward Problem.  dissertation, University of Muenster, 2019.
  * [38] G. Dassios and F. Kariotou, “On the Geselowitz formula in biomagnetics,” Quarterly of Applied Mathematics, vol. 61, no. 2, pp. 387–400, 2003.
  * [39] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, R. Kornhuber, M. Ohlberger, and O. Sander, “A generic grid interface for parallel and adaptive scientific computing. part ii: Implementation and tests in DUNE,” Computing, vol. 82, pp. 121–138, 01 2008.
  * [40] J. de Munck and M. Peters, “A fast method to compute the potential in the multisphere model (EEG application),” IEEE Transactions on Biomedical Engineering, vol. 40, no. 11, pp. 1166–1174, 1993.
  * [41] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Physics in Medicine and Biology, vol. 32, pp. 11–22, jan 1987.
  * [42] 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.
  * [43] J. C. de Munck, C. H. Wolters, and M. Clerc, EEG and MEG: forward modeling, p. 192–256.  Cambridge University Press, 2012.
  * [44] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009.
  * [45] P. Bastian, K. Birken, K. Johannsen, S. Lang, N. Neuss, H. Rentz-Reichert, and C. Wieners, “UG—a flexible software toolbox for solving partial differential equations,” Comp. Vis. Sci, vol. 1, pp. 27–40, 07 1997.
  * [46] W. Rudin, Functional Analysis.  McGraw-Hill, 2 ed., 1991.
  * [47] C. Wolters, A. Anwander, G. Berti, and U. Hartmann, “Geometry-adapted hexahedral meshes improve accuracy of finite-element-method-based EEG source analysis,” IEEE transactions on bio-medical engineering, vol. 54, pp. 1446–53, 08 2007.