# shamo: A tool for electromagnetic modeling, simulation and sensitivity
analysis of the head

Martin Grignard , Christophe Geuzaine , Christophe Phillips  M. Grignard is
with the GIGA CRC In-Vivo Imaging, University of Liège, Liège, Belgium
(e-mail: mar.grignard@uliege.be)C. Geuzaine is with the Department of
Electrical Engineering and Computer Science, University of Liège, Liège,
Belgium (e-mail: cgeuzaine@uliege.be)C. Phillips is with the GIGA CRC In-Vivo
Imaging, University of Liège, Liège, Belgium (e-mail<EMAIL_ADDRESS>0000-0001-5549-1861 0000-0001-9970-358X 0000-0002-4990-425X

(Submitted on February 14, 2022)

Accurate electromagnetic modeling of the head of a subject is of main interest
in the fields of source reconstruction and brain stimulation. Those processes
rely heavily on the quality of the model and, even though the geometry of the
tissues can be extracted from magnetic resonance images (MRI) or computed
tomography (CT), their physical properties such as the electrical conductivity
are difficult to measure with non intrusive techniques. In this paper, we
propose a tool to assess the uncertainty in the model parameters, the tissue
conductivity, as well as compute a parametric forward models for
electroencephalography (EEG) and transcranial direct current stimulation
(tDCS) current distribution.

##  1 Introduction

Accurate electromagnetic modeling of the head is of main interest for
electrophysiological source reconstruction techniques (EEG/MEG) and brain
stimulation
(tDCS/tMS). Such modeling must capture both the spatial distribution of the
tissues and their physical properties like their electrical conductivity. The
former can be extracted from different anatomical imaging techniques as
magnetic resonance imaging (MRI) and computed tomography (CT) but the latter
is very difficult to measure in vivo on a subject basis, even if some
properties can be derived from specific MRI sequences [Tuch et al., 2001, Wu
et al., 2018].

Anatomically realistic models must therefore rely on values of physical
parameters reported in the literature. The electrical conductivity of
biological tissues have been studied since the last century [Burger and
Milaan, 1943, Geddes and Baker, 1967, Gabriel et al., 1996a, b, c, Latikka et
al., 2001, Goncalves et al., 2003] and new methods are still published to
measure them accurately for each subject [Akalin Acar et al., 2016]. The
reported values have been shown to vary both inter- and intra-subject due to
temperature, health status, age, or depending on the acquisition method and
environment, e.g. in vivo vs. ex vivo [McCann et al., 2019].

Due to the variability in the published values, the influence of the chosen
parameters and geometric models on the results of the simulations have been
studied for the past decades and shown to induce erroneous electric field and
potential estimations [Haueisen et al., 1995, 1997, Vallaghe and Clerc, 2009,
Jochmann et al., 2011, Montes-Restrepo et al., 2014, Cho et al., 2015, Akalin
Acar and Makeig, 2013, Wolters et al., 2006, Vorwerk et al., 2019, Saturnino
et al., 2019]. Errors in the localisation of the reconstructed dipoles of up
to $20$ mm have been reported for basal brain locations [Lanfer et al., 2012,
Akalin Acar and Makeig, 2013]. Indeed, inaccuracies on the physical parameters
directly result in errors in the forward models, and thus in the reconstructed
sources localization or current flow.

In order to mitigate the variability in the results, different models for the
skull have been proposed [Sadleir and Argibay, 2007, Dannhauer et al., 2011,
Lanfer et al., 2012, Montes-Restrepo et al., 2014] since it acts as an
electrical insulator, in EEG and tDCS, due to its low conductivity compared to
the other tissues. While it is generally modeled as a single compartment,
partly due to the fact that further segmenting it into spongy and compact bone
is still not included in most of automated segmentation pipelines, models
differentiating these two compartments have recently been proposed [Puonti et
al., 2020, Taberna et al., 2021]. Conductivity tensors have also been
considered where the radial and tangential conductivity differ [Fuchs et al.,
2007].

The same approach applied to white matter lead to different models of its
anisotropy which have been first correlated with the water self diffusion
tensor derived from diffusion weighted MR images (DWI) [Tuch et al., 2001].
Later, the equilibrium, volume fraction and electrochemical models have been
proposed [Wu et al., 2018]. The influence of such anisotropy on EEG forward
and inverse problems have also been studied [Güllmar et al., 2010, Bashar et
al., 2010]. Conductivity tensor imaging is still an open topic with promising
advances [Ziegler et al., 2014, Sajib et al., 2016, Katoch et al., 2019].

In the past, sensitivity analysis have been mainly conducted as the final goal
of the studies. However, quantifying uncertainty in individual models could
help better understand the observed inter-subject variability in brain
stimulation and source reconstruction in broader studies. This is why we
introduce shamo, a Python open source
package111https://github.com/CyclotronResearchCentre/shamo dedicated to
stochastic electromagnetic modeling of the head and sensitivity analysis of
the results.

This toolbox aims at providing a unique solution for electromagnetic head
modeling in both source reconstruction and brain stimulation problems. While
tools already exist for each of these fields separately, for example
Brainstorm [Tadel et al., 2011] or MNE [Gramfort et al., 2013] for EEG and
SimNIBS [Thielscher et al., 2015] or ROAST [Huang et al., 2019] for tDCS to
only name a few, shamo offers an integrated solution. Moreover shamo provides
a single, easily extendable, API to perform mesh generation, simulation and
sensitivity analysis.

To highlight the mechanisms involved in our package and demonstrate its
usability and flexibility on actual cases, we apply it to the EEG forward
problem and to trans-cranial direct current stimulation (tDCS) simulation.
Both analyses are performed on a realistic finite element model (FEM) derived
from the MIDA model [Iacono et al., 2015]. To evaluate the impact of different
geometries for the skull, we build three different models, considering either
one, two or thee layers for the skull, with different electrical conductivity
values for the inner and outer tables for the latter.

The sensitivity is then assessed through the computation of Sobol indices
[Sobol, 2001]. The random input parameters considered are the values for the
electrical conductivity of the tissues. To model their probability density
functions, we use the truncated normal distribution published in the recent
review from McCann et al. [2019] as well as a unique uniform distribution, in
a worst case scenario. In the process, we compute surrogate models that, for
the EEG forward problem, results in a parametric leadfield matrix that can be
used to generate new forward models for any set of electrical conductivity
and, for the tDCS simulation, generates a model that can compute the current
density in a region of interest for the same ranges of electrical
conductivity.

##  2 Materials and methods

In order to simulate the current flow inside the brain, a mathematical model
is required. It must account for both the geometry of the tissues and their
properties. This section covers the model generation, its parametrization, and
sensitivity analysis.

###  2.1 Finite element model generation

To begin with, we focus on the geometrical aspect of the models, for which
several construction methods have been proposed [Hallez et al., 2007]: going
from a simple multi-shell sphere to a fully fledged finite element model
(FEM). Two key features of FEM are its ability to capture complex shapes and
to allow for anisotropic conductivity (in the form of a finite element field).
Pipelines have been developed to help researchers produce these models
[Windhoff et al., 2013, Nielsen et al., 2018, Huang et al., 2019, Vorwerk et
al., 2018]. As described by Huang et al. [2019], most of the available
solutions for automated segmentation rely either on Matlab, through SPM’s
“Unified Segmentation” tool [Ashburner and Friston, 2005] and its toolboxes or
on FSL and FreeSurfer [Smith et al., 2004, Fischl et al., 2004]. The ensuing
model generation step is generally tied to this specific segmentation method.
Unfortunately this prevents geometries for atypical non-healthy subjects or
simply to include other tissue types.

In shamo, we consider a FEM approach and mesh generation but eschew the
segmentation step. In effect the mesh is produced directly from a segmented
volume, i.e. where voxels are labeled as being of one of any number of tissue
classes. This allows us to work with more intricate structures and even to
model atypical cases, e.g. with prosthesis or abnormal tissue distributions
(lesions, tumor, resection,…), by using manually segmented volumes or custom
automated segmentation pipelines.

For this work, we start from the multimodal imaging-based detailed anatomical
model of the human head and neck (MIDA) [Iacono et al., 2015]: a $350\times
480\times 480$ matrix of $0.5$ $\times$ $0.5$ $\times$ $0.5$ mm3 voxels
including $116$ different structures. Based on this model, we define three
geometries with 5 to 7 different tissues (see Table 1), differing in how the
skull is modeled.

First we merge the structures of the MIDA model to keep only the main head
tissues: white matter, gray matter, cerebrospinal fluid, scalp and the
different parts of the skull. Then in model 1, the upper part of the skull is
represented as a single isotropic volume; for model 2, the upper part of the
skull is separated into spongia and compacta; for the model 3, the latter is
further divided into outer and inner tables. Note though that the lower part
of the skull is the same for all our models and is modelled as a homogeneous
tissue class. The resulting models are illustrated in Figure 1.

To generate the FEM tetrahedral mesh, we use CGAL [The CGAL Project, 2020]
with the labeled image of model 3. The resulting mesh, with $1.355\times
10^{6}$ tetrahedra, then serves as a base for the other 2 models that only
require the merging of some skull sub-compartments. This merging is performed
with Gmsh [Geuzaine and Remacle, 2009], which is also used to annotate the
mesh by specifying the names of the tissues and adding the electrodes on the
scalp. For the EEG forward problem we consider the $63$ electrodes of the
international 10-10 system [Nuwer, 2018] including the fiducial markers for
the nose, the left and right ears and the inion. The latter is considered as
the reference for the rest of this work. For tDCS, we use a subset of these
electrodes: P3, TP9, C3, P1 and O1) where P3 is the current injector.

Table 1: The tissues used in this work with the parameters of the
corresponding electrical conductivity distributions from McCann et al. [2019].
The last three columns show which tissues are included in each model. Fig. 1:
Sagittal cuts of (a) the segmented images for model 1 where skull is a single
isotropic compartment, (b) for model 2 with spongy and compact bone
differentiated, (c) for model 3 where the outer and inner tables of the skull
are differentiated and (d) the mesh corresponding to model 3. In (a), (b) and
(c), the lower part of the skull is the same.

###  2.2 Electromagnetic modeling

Since capacitive effects can be neglected in the brain tissues for the
frequency range involved in brain activity [Plonsey and Heppner, 1967], the so
called quasi-static approximation applies: the electromagnetic fields at time
$t$ only depend on the active sources at this time. In such conditions,
Maxwell’s equations reduce to a generalized Poisson equation [Malmivuo and
Plonsey, 1995, Hallez et al., 2007] that provides a relationship between the
electric potential in any point of a volume conductor and the current sources.
We first define the current density $\bm{j}$ (Am-2) and the source volume
current density $\rho_{\textrm{s}}$ (Am-3) [Schimpf, 2007]. They are linked
together by the expression

$\bm{\nabla}\cdot\bm{j}=\rho_{\textrm{s}}.$ (1)

The current density $\bm{j}$ is linearly related to the electric field
$\bm{e}$ (Vm-1) through Ohm’s law

$\bm{j}=\sigma\bm{e}$ (2)

where $\sigma$, the conductivity, can be a tensor field. Anisotropy of the
white matter has been shown to influence source reconstruction [Haueisen et
al., 2002, Güllmar et al., 2010] and the pipeline allows for anisotropic
tissues yet, in this work, $\sigma$ is considered isotropic because no
diffusion weighted images (DWI) is available in the MIDA data. The quasi-
static conditions described above allow us to write the relationship between
the electric field and the electric potential field $v$ (V) as

$\bm{e}=-\bm{\nabla}v.$ (3)

Then combining Equations (1), (2) and (3) leads to the generalized Poisson
equation:

$\bm{\nabla}\cdot(\sigma\bm{\nabla}(v))=-\rho_{\textrm{s}}.$ (4)

Finally a homogeneous Neumann boundary condition is set on the interface
between the conductor volume and the air, and a Dirichlet condition is added
to set the reference electrode. We use GetDP [Geuzaine, 2007] as a solver for
the finite element model described in section 2.1. GetDP is an industrial
grade solver and is at the heart of some popular tools for head
electromagnetic modeling [Huang et al., 2019] and our implementation of the
forward problem have already been analytically validated by Ziegler et al.
[2014].

###  2.3 Electrical conductivity of tissues

As stated in Equation (4), electrical conductivity plays a major role in the
computation of the electric potential and the other related fields.
Unfortunately, determining the exact electrical conductivity $\sigma$ (Sm-1)
of the biological tissues in a non-intrusive manner is an open issue. Multiple
methods have been developed to measure it either in vitro, ex vivo or even in
vivo but struggle to provide an accurate and reliable value [Burger and
Milaan, 1943, Geddes and Baker, 1967, Gabriel, 1996, Gabriel et al., 1996a, b,
c, Latikka et al., 2001, Goncalves et al., 2003]. McCann et al. [2019]
reviewed the values acquired with different techniques and under specific
conditions and derived a realistic underlying probability distribution in the
form of a truncated normal distribution for the main tissues composing the
head. We use these distributions with the 3 models described in Section 2.1,
which we label Models 1a, 2a and 3a.

In addition, we define a uniform distribution spanning the whole range of the
reported conductivity values for all the tissue classes, and refer to it as
the extended distribution (EXT). This represents a worst case scenario with no
prior information on the conductivity of the different tissues. This uniform
distribution is also used with the 3 FEMs, which we label Models 1b, 2b and
3b. The range and distribution of conductivity values for the tissues
considered in each model are summarized in Table 1.

The goal of the sensitivity analysis, see Section 2.5, is to determine the
parameters that drive the variability in the results. Comparing the
sensitivity to the two sets of conductivity values, realistic and extended,
will allow us to assess the effect of the prior knowledge added by the
truncated normal distributions, especially for the tissues with the narrowest
distributions.

For the sake of clarity, we use lowercase letters to indicate a known
deterministic value whereas uppercase letters refer to random values. For
instance, the tissues conductivities are denoted by the vector
$\bm{\sigma}=[\sigma_{1},\>\sigma_{2},\>\dots,\>\sigma_{d}]$ and
$\bm{\Sigma}=[\Sigma_{1},\>\Sigma_{2},\>\dots,\>\Sigma_{d}]$ corresponds to
the vector containing the random parameters modeled by the distributions. Thus
for each geometry $i$, we consider two sets of conductivity distributions:
either those introduced in [McCann et al., 2019], giving
$\bm{\Sigma}_{\textrm{a}}=[\Sigma_{\textrm{WM}},\>\Sigma_{\textrm{GM}},\>\dots,\>\Sigma_{\textrm{SCP}}]$,
or the same extended distribution for each of the tissues, i.e.
$\bm{\Sigma}_{\textrm{b}}=[\Sigma_{\textrm{EXT}},\>\Sigma_{\textrm{EXT}},\>\dots,\>\Sigma_{\textrm{EXT}}]$
as shown in Table 1.

###  2.4 EEG and tDCS forward problem

When carrying out EEG source reconstruction analysis, one attempts to recover
the underlying brain activity inducing the observed signal at the scalp level.
This electrical activity is generally modeled by one or more equivalent
current dipoles characterized by their coordinates in space
$\bm{r}=[r_{x},\>r_{y},\>r_{z}]$ and their dipole moment
$\bm{p}=[p_{x},\>p_{y},\>p_{z}]$ (Am). In practice, a set of discrete sources
is considered rather than the full continuous volume of the gray matter. This
set is called a source space and defines potential dipole locations.

The relation between the source space containing $n$ sources and the electric
potential measured on $m$ electrodes at the scalp level is given by the
expression

$[l]\cdot\bm{s}+\bm{\varepsilon}=\bm{v},$ (5)

where
$\bm{s}=[p_{1}^{(x)},\>p_{1}^{(y)},\>p_{1}^{(z)},\>\dots,\>p_{n}^{(x)},\>p_{n}^{(y)},\>p_{n}^{(z)}]^{\mathsf{T}}$
is the source vector and $p_{j}^{(k)}$ is the dipole moment of the source
located in the $j$-th site along $k$-axis,
$\bm{\varepsilon}=[\varepsilon_{1},\>\dots,\>\varepsilon_{m}]^{\mathsf{T}}$
and $\bm{v}=[v_{1},\>\dots,\>v_{m}]^{\mathsf{T}}$ are respectively the
additive noise component vector and the vector of electrodes potentials (V),
and $[l]$ is equally referred to as the ”leadfield” or gain matrix. This
matrix looks like this

$[l]=\begin{bmatrix}l_{1,1}^{(x)}&l_{1,1}^{(y)}&l_{1,1}^{(z)}&\dots&l_{1,n}^{(x)}&l_{1,n}^{(y)}&l_{1,n}^{(z)}\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\\
l_{m,1}^{(x)}&l_{m,1}^{(y)}&l_{m,1}^{(z)}&\dots&l_{m,n}^{(x)}&l_{m,n}^{(y)}&l_{m,n}^{(z)}\end{bmatrix},$
(6)

where each element $l_{i,j}^{(k)}$ corresponds to the electric potential $v$
measured on the $i$-th electrode due to a current dipole with unitary dipole
moment located on the $j$-th site and oriented along $k$-axis (VA-1m-1). This
model encompasses all the geometric information and the physical properties of
the head tissues. In the rest of this paper, the notation $[l(\bm{\sigma})]$
is used when highlighting the dependencies of the leadfield matrix on the
values set for the electrical conductivity.

Following the method described by Weinstein et al. [2000] based on the
reciprocity principle, we actually have to solve the tDCS forward problem in
order to generate the EEG leadfield. Indeed, to compute $[l]$ on an element
basis, the reciprocity principle states that to estimate the voltage
difference between two points due to a single current dipole, one needs to
compute the electric field $\bm{e}$ at the coordinates of the dipole resulting
from the injection of a $1$ A current $i$ between the two points, which is the
definition of a tDCS forward problem.

$v_{1}-v_{2}=\frac{\bm{e}\cdot\bm{p}}{i}$ (7)

The technique then consists in setting a reference electrode, iteratively
injecting a $1$ A current through the $m$ active electrodes, and estimating
the electric field on the source space in the $i$-th row of the matrix $[l]$.
This step of the process is achieved in GetDP with the generalized minimal
residual method (GMRES) configured with a tolerance of $10^{-8}$ and an
incomplete factorization (ILU) preconditioner.

Given that the current sources should only exist in the gray matter, that the
mesh for that tissue is made of about $368000$ tetrahedra, and that we have a
setup of $59$ active electrodes, the whole leadfield matrix
$[l]_{\textrm{full}}$ would theoretically have a size of $59\times(3\times
368000)$, which is too large for practical use. Therefore we arbitrarily fix
the average interval between two sources at $7.5$ mm, resulting in $2127$
sources and a leadfield matrix $[l]\in\mathbb{R}^{59\times 6381}$ which is
more manageable. This source-to-source distance influences both the the
computational resources required to perform the source reconstruction, since
it is directly linked to the size of the leadfield matrix, and the number of
potentially reconstructed dipoles. In their paper, Michel and Brunet [2019]
state that the definition of the spatial resolution is a sensitive problem but
that increasing it does not lead to a linear increase of the accuracy. In
fact, the accuracy has a limit due to the fact that the amount of information
provided to the inverse problem remains constant since the number of
electrodes is fixed.

###  2.5 Sensitivity analysis

As defined by Saltelli [2008], sensitivity analysis is the study of how
variation in the input parameters of a process influences the variation in the
output. In this field, two cases are differentiated. The local sensitivity
focuses on the uncertainty at a specific coordinate of the parameters space
$\Omega$ whereas the global sensitivity captures the variation across the
whole space.

One of the most used and studied global sensitivity analysis techniques is the
computation of the so called Sobol indices [Sobol, 2001, Saltelli et al.,
2010]. Let us consider a model $Y=Y(\bm{X})$ where $Y$ is the random output
variable and $\bm{X}=[X_{1},\dots,X_{n_{\textrm{p}}}]$ is the vector of
$n_{\textrm{p}}$ random input variables. The first and total order Sobol
indices for the $i$-th input variable, $s_{i}$ and $s_{i}^{\textrm{(t)}}$ are
defined by

$\displaystyle s_{i}$
$\displaystyle=\frac{\mathbb{V}_{X_{i}}(\mathbb{E}_{\bm{X}_{\setminus
i}}(Y\;|\;X_{i}))}{\mathbb{V}(Y)},$ (8) $\displaystyle s_{i}^{\textrm{(t)}}$
$\displaystyle=\frac{\mathbb{E}_{\bm{X}\setminus
i}(\mathbb{V}_{X_{i}}(Y\;|\;\bm{X}_{\setminus i}))}{\mathbb{V}(Y)},$ (9)

where $\bm{X}_{\setminus i}$ is the vector of all the random inputs but
$X_{i}$, $s_{i}$ corresponds to the variance in the output explained by
$X_{i}$ alone, $\mathbb{V}_{X_{i}}(\mathbb{E}_{\bm{X}_{\setminus
i}}(Y\;|\;X_{i}))$ is the variance explained by the $i$-th parameter, also
referred to as its main effect, and $s_{i}^{\textrm{(t)}}$ is the output
variance explained by $X_{i}$ and all its interactions with the other input
parameters.

To compute Sobol indices, we follow the method presented by Saltelli et al.
[2010] and implemented in the python package SALib [Herman and Usher, 2017]
that provides a way to compute both $s_{i}$ and $s_{i}^{(t)}$ from the same
set of evaluations of the model, thus reducing the amount of computation
required. This technique relies on $n_{d}$ observations
$\\{(y_{i}^{\textrm{(d)}},\bm{x}_{i}^{\textrm{(d)}},\;i=1,\dots,n_{d}\\}$
where each $y_{i}^{\textrm{(d)}}=y(\bm{x}_{i}^{\textrm{(d)}})$ is the output
of the model for a set of inputs
$\bm{x}_{i}^{\textrm{(d)}}=[x_{i,1}^{\textrm{(d)}},\dots,x_{i,n_{p}}^{\textrm{(d)}}]$.
Let us define
$\bm{y}^{\textrm{(d)}}=[y_{1}^{\textrm{(d)}},\dots,y_{n_{d}}^{\textrm{(d)}}]^{\mathsf{T}}$
the vector of outputs and
$[x]^{\textrm{(d)}}=[\bm{x}_{1}^{\textrm{(d)}},\dots,\bm{x}_{n_{d}}^{\textrm{(d)}}]^{\mathsf{T}}$
the matrix of inputs.

The matrix $[x]^{\textrm{(d)}}$ is built of $n_{\textrm{p}}+2$ sub-matrices:
$[a]$, $[b]$ and the matrices $[a_{b}]^{(i)}$ where all the columns are the
same as in $[a]$ except the $i$-th one coming from $[b]$. All these matrices
have $n_{r}$ rows and $n_{\textrm{p}}$ columns. The input vectors
$\bm{x}_{i}^{\textrm{(d)}}$ composing the independent matrices $[a]$ and $[b]$
are drawn from the parameters space $\Omega$ using the Saltelli extension of
Sobol quasi-random sequence [Sobol, 1967, 1976]. Such sequences are described
in section 2.6.

Based on these samples, the numerators of Equations (8) and (9) are computed
with

$\begin{split}&\mathbb{V}_{X_{i}}(\mathbb{E}_{\bm{X}_{\setminus
i}}(Y\;|\;X_{i}))\\\
&\qquad=\frac{1}{n_{r}}\sum_{j=1}^{n_{r}}\bm{y}([b])_{j}\left(\bm{y}([a_{b}])_{j}^{(i)}-\bm{y}([a])_{j}\right),\end{split}$
(10)

and

$\begin{split}&\mathbb{E}_{\bm{X}\setminus
i}(\mathbb{V}_{X_{i}}(Y\;|\;\bm{X}_{\setminus i}))\\\
&\qquad=\frac{1}{2n_{r}}\sum_{j=1}^{n_{r}}\left(\bm{y}([a])_{j}-\bm{y}([a_{b}])_{j}^{(i)}\right)^{2}.\end{split}$
(11)

###  2.6 Surrogate model

The computation of the sensitivity indices described in section 2.5 requires a
large number of model evaluations. When the estimation of the actual model
(here the computation of the leadfield matrix) is computationally heavy, a
simpler model, referred to as the surrogate model, can be used instead. This
simpler version must behave almost like if it were the real one but its
evaluation should require less computing power.

Building such a model is the goal of all the supervised learning techniques.
Those methods start from a set of $n_{d}$ observations
$\bm{y}^{\textrm{(d)}}=[y_{1}^{\textrm{(d)}},\>\dots,\>y_{n_{d}}^{\textrm{(d)}}]^{\mathsf{T}}$
of the actual model at different points of the parameters space
$[x]^{\textrm{(d)}}=[\bm{x}_{1}^{\textrm{(d)}},\>\dots,\>\bm{x}_{n_{d}}^{\textrm{(d)}}]^{\mathsf{T}}$
where $y_{i}^{\textrm{(d)}}=y(\bm{x}_{i}^{\textrm{(d)}})$ with $y(\bm{x})$ the
real model. From this relatively small amount of evaluations of the model, the
surrogate model $\hat{y}(\bm{x})$ is built so that
$\hat{y}=\hat{y}(\bm{x})\approx y(\bm{x})$ for any vector $\bm{x}\in\Omega$
that is not in the training set $[x]^{\textrm{(d)}}$.

The first step for building the surrogate model is then to draw $n_{d}$
vectors $\bm{x}_{i}^{\textrm{(d)}}$ to build the matrix $[x]^{\textrm{(d)}}$.
This can be performed with various methods but here we consider quasi-random
sequences. Those sequences, compared to real random sequences, take into
account the previous points that have been drawn. They are used to cover the
space as efficiently as possible. In section 2.5 the Saltelli extension of
Sobol sequence is used to define the coordinates and here, to produce the
training set for the surrogate model, we use a Halton sequence [Halton, 1960]
as implemented in chaospy [Feinberg and Langtangen, 2015].

In shamo, the generation of the surrogate model is carried out with “Gaussian
Processes Regression” (GPR) [Rasmussen and Williams, 2006]. Let us define the
notation for a multivariate normal distribution
$\mathcal{N}(\bm{\mu},[\gamma])$ where $\bm{\mu}=[\mu_{1},\dots,\mu_{n_{p}}]$
is the vector of means along each axis and $[\gamma]$ is the covariance matrix
where each $\gamma_{i,i}$ is the variance of the $i$-th random parameter and
the elements $\gamma_{i,j}$ are the correlation between the $i$-th and the
$j$-th variables.

To predict the $n_{t}$ values
$\bm{y}^{(t)}=[y_{1}^{(t)},\dots,y_{n_{t}}^{(t)}]^{\mathsf{T}}$ on the test
points $[x]^{(t)}$, GPR handles the problem as Bayesian inference. Under these
conditions, the learning samples are treated as random variables following a
multivariate normal distribution
$P(\bm{y}^{\textrm{(d)}}\;|\;[x]^{\textrm{(d)}}=\mathcal{N}(\mu^{\textrm{(d)}},[\gamma]^{\textrm{(d)}})$.
Here, the mean of this distribution is set to the mean of the learning
outputs. To consider the test points, this expression becomes

$\begin{split}P&(\bm{y}^{\textrm{(d)}},\bm{y}^{(t)}\;|\;[x]^{\textrm{(d)}})\\\
&=\mathcal{N}\left(\mu^{\textrm{(d)}},\begin{bmatrix}[\gamma]^{(t)}&[\gamma]^{(t,d)}\\\
[\gamma]^{(d,t)}&[\gamma]^{\textrm{(d)}}+\epsilon[i]\end{bmatrix}\right),\end{split}$
(12)

with $\epsilon$ an added noise.

Next, the conditional distribution
$P(\bm{y}^{(t)}\;|\;\bm{y}^{\textrm{(d)}},[x]^{\textrm{(d)}})=\mathcal{N}(\mu^{*},[\gamma]^{*})$
is obtained with

$\displaystyle\mu^{*}$
$\displaystyle=\mu^{\textrm{(d)}}+[\gamma]^{(t,d)}([\gamma]^{\textrm{(d)}})^{-1}(\bm{y}^{\textrm{(d)}}-\mu^{\textrm{(d)}}),$
(13) $\displaystyle[\gamma]^{*}$
$\displaystyle=[\gamma]^{(t)}-[\gamma]^{(t,d)}([\gamma]^{\textrm{(d)}})^{-1}[\gamma]^{(d,t)}.$
(14)

Finally, the mean values $\mu_{i}^{*}$ and the standard deviation
$\gamma_{i}^{*}=[\gamma]_{i,i}^{*}$ are obtained by the marginalisation of
each random variable. The values $\mu_{i}^{*}$ are the predictors
corresponding to the test points.

During the training step, the hyper-parameters of the kernel are optimised by
maximising the log-marginal likelihood (LML) [Schirru et al., 2011]. When the
model outputs more than one scalar, the process can be applied separately to
each of the outputs, giving one Gaussian process by output variable. Here, we
use the implementation of the GPR from scikit-learn [Pedregosa et al., 2011]
and follow the recommendations from Chen et al. [2016]. Thus, the regression
part of the GPR is set to the mean of the output variable and the kernel is
obtained by the product of a constant kernel and a stationary Matérn kernel
with the smoothness parameter $\nu=2.5$, thus resulting in the covariance
function

$\begin{split}k(\bm{x}_{1},\bm{x}_{2})=&\left(1+\frac{\sqrt{5}}{l}d(\bm{x}_{1},\bm{x}_{2})+\frac{5}{3l}d(\bm{x}_{1},\bm{x}_{2})r\right)\\\
&\cdot\textrm{exp}\left(-\frac{\sqrt{5}}{l}d(\bm{x}_{1},\bm{x}_{2})\right).\end{split}$
(15)

##  3 Applications

We demonstrate the application of shamo on EEG and tDCS forward problems.

###  3.1 EEG forward problem

As described in section 2.4, the computation of the EEG forward model is of
main interest in source reconstruction but is highly dependent on the geometry
and the physical properties of the tissues.

To build the surrogate model, we generate a set of leadfield matrices
$[l(\bm{\sigma^{(i)})}]$ for $100$ conductivity vectors $\bm{\sigma}^{(i)}$
drawn from the parameters space $\Omega$ using a Halton sequence. This step
results in a leadfield matrix where each element is actually a Gaussian
process, which gives us the ability to quickly construct any new matrix
$[\hat{l}(\bm{\sigma})]$ for a specific conductivity vector $\bm{\sigma}$.

The sensitivity indices defined in Equations (8) and (9) are only valid for a
model with a single scalar output. Therefore we choose to study the
sensitivity of the whole matrix to the values of $\bm{\sigma}$ with a distance
measure $m(\bm{\sigma})$ relative to a reference leadfield matrix
$[l]_{\textrm{ref}}$, obtained with a fixed
$\bm{\sigma}=\bm{\sigma}_{\textrm{ref}}$

$m(\bm{\sigma})=\left\lVert[l(\bm{\sigma})]-[l]_{\textrm{ref}}\right\rVert_{\mathsf{F}},$
(16)

where $\bm{\sigma}_{\textrm{ref}}$ is the mean value for each tissue (See
Table 1) and and $\left\lVert\dots\right\rVert_{\mathsf{F}}$ is the Frobenius
norm.

A surrogate model $\hat{m}(\bm{\sigma})$ is thus built for this function based
on the same training data as the parametric matrix introduced above. Next, the
first and total order Sobol indices are computed from two sets of $40000$
evaluations of the Gaussian process for the six models of this study, as
defined in Section 2.3. The resulting indices are shown in Figure 2.

Clearly, for both the truncated distributions and the extended uniform ones,
the parameter with the largest influence on the metric is the gray matter
conductivity $\sigma_{\textrm{GM}}$. Whether one uses the narrow truncated
normal (models 2a/3a) or the extended uniform (models 2b/3b) distributions for
the compact skull and the outer and inner tables has little influence on the
Sobol indices.

Another interesting point is the increasing influence of the CSF conductivity
when the uniform distribution is used. While both its first and total order
Sobol indices in models 1a to 3a are very small, the value of
$s_{\textrm{CSF}}^{(t)}$ for models 1b to 3b are non negligible meaning there
are interactions between parameters.

Fig. 2: First (%) and total order Sobol indices of the metric $m(\bm{\sigma})$
for each tissue of each model. In the left column, the values for
$\bm{\Sigma}$ are the truncated normal distributions from McCann et al. [2019]
and int the right column, the extended uniform distribution is used.

###  3.2 Evaluation of the electrodes potential

To illustrate the actual effect of the sensitivity described in the above
application, we calculate the electrical potential measured on the scalp due
to a single left frontal source located $17.8$mm under F3. Model 3 with the
$\bm{\sigma}_{\textrm{ref}}$ was used as a reference (Figure 3b) then the
conductivity of GM was also set to the minimal and maximal value found in the
literature (See Table 1) leading to slightly modified electrical potential
scalp maps (Figure 3a, c). The scalp map differences of the latter two with
the reference one is shown in Figure 3 d, e.

Fig. 3: (a) the scalp potential (V) computed for
$\sigma_{\textrm{GM}}=\sigma_{\textrm{GM,min}}=0.06$ (Sm-1), (b) the scalp
potential (V) obtained for
$\sigma_{\textrm{GM}}=\sigma_{\textrm{GM,ref}}=0.466$ (Sm-1) and (c) the scalp
potential (V) measured with
$\sigma_{\textrm{GM}}=\sigma_{\textrm{GM,max}}=2.47$. (d) and (e) show the
difference between the computed scalp potential and the reference one (V)
respectively for $\sigma_{\textrm{GM}}=\sigma_{\textrm{GM,min}}$ and
$\sigma_{\textrm{GM}}=\sigma_{\textrm{GM,max}}$.

###  3.3 Transcranial direct current stimulation (tDCS) simulation

Using the same formulation as for the EEG forward problem, we can model tDCS
and obtain the current density, electric potential and electric field accross
the whole head.

To illustrate this aspect, we consider a HD-tDCS experiment where electrode P3
was set as a 4 mA injector and electrodes TP9, C3, P1 and O1 were set to
ground. We used the mesh from model 3 and the truncated normal conductivity
distributions as in model 3a for this simulation. As a metric to assess the
sensitivity of the model with regard to the conductivity values, we chose the
mean of current density norm in a small region of 368 mm3 located $22.6$ mm
under CP3. The results of these simulations are shown in Figure 4. As an extra
feature for researchers in neuroscience, the estimated fields can also be
directly exported as a standard multidimensional NifTI image.

Fig. 4: (a) A view of the scalp potential induced by the set electrodes
surrounded in white, (b) the first (%) and total order Sobol indices of the
mean current density in the ROI for each tissue conductivity, (c) a cut of the
current density inside the head resulting from the injection of current inside
the reference model where $\bm{\sigma}=\bm{\sigma}_{\textrm{ref}}$ and (d) the
same information in the form of a NIFTI file with a representation of the
small mask used to compute the mean of the norm of the current density.

As visible in Figure 4, current density is highest in the scalp tissue between
the electrodes but also spreads diffusely throughout the head volume. The
tissues with the highest Sobol index are the scalp, followed by the spongy
compartment of the skull. This comes as no surprise as electrical current
follows the path of least resistance.

##  4 Discussion

The pipeline presented in this work uses several well established methods.
Here, we discuss the added value of such tool and technical choices and
compare it to established solutions such as SimNIBS [Thielscher et al., 2015]
or BrainStorm [Tadel et al., 2011].

First, the generation of a realistic subject specific head model is generally
a tedious task involving the segmentation of the head volume, based on MR or
CT images, then the delineation of the tissue interfaces. This step often
requires further ad hoc cleaning to ensure the surfaces are two dimensional
manifold, i.e. they are completely closed and thus have an inside and an
outside. Then, those surfaces must be integrated into a single mesh and the so
defined volumes filled with tetrahedral elements.

Several automated pipelines are now available is the most common toolboxes
[Thielscher et al., 2015, Huang et al., 2019, Nielsen et al., 2018] using
either SPM [Ashburner and Friston, 2005] or FSL/FreeSurfer [Smith et al.,
2004, Fischl et al., 2004] to segment the structural images. While this
greatly simplifies model generation, since it takes care of both the
segmentation and meshing steps, it also prevents the user from using non
standard segmentations.

Here the FEM mesh is directly built from a 3D image made with labeled voxels.
This approach not only provides more control on the refinement of the final
mesh but also allows the any number of tissues, even with complex
configurations. In practice, this opens the path for more specific cases such
as modeling prosthetic, metallic plates or any other unusual head geometry.
Moreover one can freely choose the segmentation tool to use or could even
proceed manually for difficult cases.

As explained previously, the model presented in this paper does not consider
the anisotropy of white matter because no tensor information is provided with
the original data. However any kind of field can be included in the model and
handled by the solver, GetDP, without extra burden from the point of view of
the user. Besides GetDP already powers other similar tools like SimNIBS
[Thielscher et al., 2015] and ROAST [Huang et al., 2019] and the forward
problem implementation used in shamo have been validated by [Ziegler et al.,
2014].

Here we only considered two conductivity dependent electrostatic problems, EEG
and tDCS. Nevertheless, thanks to GetDP’s simple formalism, shamo could be
directly extended to electromagnetic modeling processes, including MEG and
TMS. Indeed this only requires the definition of the equations related to the
problem, which are explicitly stored in a problem file for GetDP. Nowadays
other tools, like SimNIBS, also include modules for the quantification of the
effect of tissue conductivity uncertainty [Saturnino et al., 2019], showing
the importance of such analysis for tDCS applications. Still, beyond
electrical conductivity, the sensitivity analysis could also focus on other
parameters, e.g. the injected current in HD-tDCS, through python classes
available in shamo to expose the needed parameters. Importantly, by
definition, the Sobol indices rely on a single scalar output function, whose
choice is thus critical but also lets the user focus on any feature of
interest. Thereby shamo’s flexible implementation allows one to define his own
processes and sensitivity analysis.

Regarding the surrogate models, we decided to use Gaussian process because
they provide information on the confidence over the solution through the
standard deviation on each predicted point. This can be used to obtain more
in-depth understanding of the model. Such regressor also has the advantage of
not requiring huge amounts of training data. Considering the fact that an
evaluation of the actual model can take several hours, reducing the number of
observations can drastically lower the computation time required. To further
reduce this time, the tool provides an easy way to evaluate each sample point
separately, thus allowing the use of high-performance computing (HPC)
equipment like computer clusters. In the present work, observations were
computed by batches of $100$, each on a single core on the CÉCI
clusters222http://www.ceci-hpc.be/.

Overall the goal of shamo is to provide a single tool to perform three major
steps, namely FEM creation, model estimation, and sensitivity analysis. This
operated with few dependencies that are all established, in an open source
software working out of the box on any major operating system or on HPC
platforms.

##  5 Conclusion

In this paper, we presented a python pipeline for accurate electromagnetic
modeling of the head which allows for sensitivity analysis and surrogate model
building, bringing together similar features as some established software,
such as SimNIBS [Thielscher et al., 2015] or ROAST [Huang et al., 2019] for
tDCS and Brainstorm [Tadel et al., 2011] or MNE [Gramfort et al., 2013] for
EEG, unified with a single API. This tool, called shamo [Grignard, 2021a], and
the full documentation [Grignard, 2021b] are available on Github under GPL-v3
license. A set of examples are also available in the form of jupyter
notebooks.

We showed a use-case for the EEG forward problem where a parametric leadfield
matrix is computed and can then be used to generate any new matrix for a
specific set of tissue conductivity values and another application to tDCS
where the current density in a certain region is obtained and can be studied
with regard to the electrical sensitivity. Those are only two possible
applications but shamo could easily be extended to magnetic stimulation or
TMS.

Considering the abstraction level of the tool and the outcome that can be
obtained from it, one can use our tool to perform finite element analysis and
sensitivity analysis without having to dig into those fields, letting the user
employ the toolset of his choice for further analysis. shamo could be used in
various studies to assess the sensitivity of the results to some parameters or
to build parametric models for complex physical fields that, otherwise, should
be evaluated every time a new value is needed.

## Data availability

The data that support the findings of this study are available from the IT’IS
foundation333https://itis.swiss/virtual-population/regional-human-models/mida-
model/ but restrictions apply to the availability of these data, which were
used under
licence444https://itis.swiss/assets/Downloads/VirtualPopulation/License_Agreements/LicenseAgreementMIDA.pdf
for the current study, and so are not publicly available. Data are however
available from the authors upon reasonable request and with permission of the
IT’IS foundation.

## Information sharing

The source code of shamo is available on
Github555https://github.com/CyclotronResearchCentre/shamo under GPLv3 license
and is fully
documented666https://cyclotronresearchcentre.github.io/shamo/index.html. It
can be installed from PyPI777https://pypi.org/project/shamo/. In addition,
jupyter notebooks are also available on
Github888https://github.com/CyclotronResearchCentre/shamo-tutorials and show
how to conduct similar studies.

## Acknowledgements

MG and CP are supported by the Fonds de la Recherche Scientifique de Belgique
(F.R.S.-FNRS), the former under Grant No. EOS 30446199, Belgium.

Computational resources have been provided by the Consortium des Équipements
de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de
Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region,
Belgium.

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * Akalin Acar and Makeig [2013] Zeynep Akalin Acar and Scott Makeig.  Effects of forward model errors on EEG source localization.  _Brain Topography_ , 26(3):378–396, 2013.  ISSN 08960267.  doi: 10.1007/s10548-012-0274-6.  URL http://link.springer.com/10.1007/s10548-012-0274-6.
  * Akalin Acar et al. [2016] Zeynep Akalin Acar, Can E. Acar, and Scott Makeig.  Simultaneous head tissue conductivity and EEG source location estimation.  _NeuroImage_ , 124:168–180, January 2016.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2015.08.032.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811915007442.
  * Ashburner and Friston [2005] John Ashburner and Karl J. Friston.  Unified segmentation.  _NeuroImage_ , 26(3):839–851, July 2005.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2005.02.018.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811905001102.
  * Bashar et al. [2010] M. R. Bashar, Y. Li, and P. Wen.  Uncertainty and sensitivity analysis for anisotropic inhomogeneous head tissue conductivity in human head modelling.  _Australasian Physical & Engineering Sciences in Medicine_, 33(2):145–152, June 2010.  ISSN 0158-9938, 1879-5447.  doi: 10.1007/s13246-010-0015-7.  URL http://link.springer.com/10.1007/s13246-010-0015-7.
  * Burger and Milaan [1943] H. C. Burger and J. B. van Milaan.  Measurements of the specific Resistance of the human Body to direct Current.  _Acta Medica Scandinavica_ , 114(6):584–607, 1943\.  ISSN 0954-6820.  doi: 10.1111/j.0954-6820.1943.tb11253.x.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/j.0954-6820.1943.tb11253.x.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.0954-6820.1943.tb11253.x.
  * Chen et al. [2016] Hao Chen, Jason L. Loeppky, Jerome Sacks, and William J. Welch.  Analysis Methods for Computer Experiments: How to Assess and What Counts?  _Statistical Science_ , 31(1):40–60, February 2016.  ISSN 0883-4237.  doi: 10.1214/15-STS531.  URL http://projecteuclid.org/euclid.ss/1455115913.
  * Cho et al. [2015] Jae-Hyun Cho, Johannes Vorwerk, Carsten H. Wolters, and Thomas R. Knösche.  Influence of the head model on EEG and MEG source connectivity analyses.  _NeuroImage_ , 110:60–77, April 2015.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2015.01.043.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811915000683.
  * Dannhauer et al. [2011] Moritz Dannhauer, Benjamin Lanfer, Carsten H. Wolters, and Thomas R. Knösche.  Modeling of the human skull in EEG source analysis.  _Human Brain Mapping_ , 32(9):1383–1399, September 2011.  ISSN 10659471.  doi: 10.1002/hbm.21114.  URL http://doi.wiley.com/10.1002/hbm.21114.
  * Feinberg and Langtangen [2015] Jonathan Feinberg and Hans Petter Langtangen.  Chaospy: An open source tool for designing methods of uncertainty quantification.  _Journal of Computational Science_ , 11:46–57, November 2015.  ISSN 18777503.  doi: 10.1016/j.jocs.2015.08.008.  URL https://linkinghub.elsevier.com/retrieve/pii/S1877750315300119.
  * Fischl et al. [2004] Bruce Fischl, David H. Salat, André J. W. van der Kouwe, Nikos Makris, Florent Ségonne, Brian T. Quinn, and Anders M. Dale.  Sequence-independent segmentation of magnetic resonance images.  _NeuroImage_ , 23 Suppl 1:S69–84, 2004.  ISSN 1053-8119.  doi: 10.1016/j.neuroimage.2004.07.016.
  * Fuchs et al. [2007] Manfred Fuchs, Michael Wagner, and Joern Kastner.  Development of volume conductor and source models to localize epileptic foci.  _Journal of Clinical Neurophysiology: Official Publication of the American Electroencephalographic Society_ , 24(2):101–119, April 2007.  ISSN 0736-0258.  doi: 10.1097/WNP.0b013e318038fb3e.
  * Gabriel et al. [1996a] C Gabriel, S Gabriel, and E Corthout.  The dielectric properties of biological tissues: I. Literature survey.  _Physics in Medicine and Biology_ , 41(11):2231–2249, November 1996a.  ISSN 0031-9155, 1361-6560.  doi: 10.1088/0031-9155/41/11/001.  URL https://iopscience.iop.org/article/10.1088/0031-9155/41/11/001.
  * Gabriel [1996] Camelia Gabriel.  Compilation of the Dielectric Properties of Body Tissues at RF and Microwave Frequencies.:.  Technical report, Defense Technical Information Center, Fort Belvoir, VA, January 1996.  URL http://www.dtic.mil/docs/citations/ADA303903.
  * Gabriel et al. [1996b] S Gabriel, R W Lau, and C Gabriel.  The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz.  _Physics in Medicine and Biology_ , 41(11):2251–2269, November 1996b.  ISSN 0031-9155, 1361-6560.  doi: 10.1088/0031-9155/41/11/002.  URL https://iopscience.iop.org/article/10.1088/0031-9155/41/11/002.
  * Gabriel et al. [1996c] S Gabriel, R W Lau, and C Gabriel.  The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues.  _Physics in Medicine and Biology_ , 41(11):2271–2293, November 1996c.  ISSN 0031-9155, 1361-6560.  doi: 10.1088/0031-9155/41/11/003.  URL https://iopscience.iop.org/article/10.1088/0031-9155/41/11/003.
  * Geddes and Baker [1967] L. A. Geddes and L. E. Baker.  The specific resistance of biological material–a compendium of data for the biomedical engineer and physiologist.  _Medical & Biological Engineering_, 5(3):271–293, May 1967.  ISSN 0025-696X.  doi: 10.1007/BF02474537.
  * Geuzaine [2007] C. Geuzaine.  GetDP: a general finite-element solver for the de Rham complex.  _PAMM Volume 7 Issue 1. Special Issue: Sixth International Congress on Industrial Applied Mathematics (ICIAM07) and GAMM Annual Meeting, Zürich 2007_ , 7(1):1010603–1010604, 2007.  doi: 10.1002/pamm.200700750.
  * Geuzaine and Remacle [2009] Christophe Geuzaine and Jean-François Remacle.  Gmsh: a three-dimensional finite element mesh generator with built-in pre- ad post-processing facilities.  _International Journal for Numerical Methods in Engineering_ , 79(11):1309–1331, 2009.  ISSN 00295981.  doi: 10.1002/nme.2579.
  * Goncalves et al. [2003] S.I. Goncalves, J.C. de Munck, J.P.A. Verbunt, F. Bijma, R.M. Heethaar, and F. Lopes da Silva.  In vivo measurement of the brain and skull resistivities using an EIT-based method and realistic models for the head.  _IEEE Transactions on Biomedical Engineering_ , 50(6):754–767, June 2003.  ISSN 1558-2531.  doi: 10.1109/TBME.2003.812164.  Conference Name: IEEE Transactions on Biomedical Engineering.
  * Gramfort et al. [2013] Alexandre Gramfort, Martin Luessi, Eric Larson, Denis A. Engemann, Daniel Strohmeier, Christian Brodbeck, Roman Goj, Mainak Jas, Teon Brooks, Lauri Parkkonen, and Matti S. Hämäläinen.  MEG and EEG data analysis with MNE-Python.  _Frontiers in Neuroscience_ , 7(267):1–13, 2013\.  doi: 10.3389/fnins.2013.00267.
  * Grignard [2021a] Martin Grignard.  shamo.  https://doi.org/10.5281/zenodo.4420811, January 2021a.  URL https://doi.org/10.5281/zenodo.4420811.
  * Grignard [2021b] Martin Grignard.  shamo documentation.  https://cyclotronresearchcentre.github.io/shamo/index.html, January 2021b.  URL https://cyclotronresearchcentre.github.io/shamo/index.html.
  * Güllmar et al. [2010] Daniel Güllmar, Jens Haueisen, and Jürgen R. Reichenbach.  Influence of anisotropic electrical conductivity in white matter tissue on the EEG/MEG forward and inverse solution. A high-resolution whole head simulation study.  _NeuroImage_ , 51(1):145–163, May 2010.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2010.02.014.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811910001825.
  * Hallez et al. [2007] Hans Hallez, Bart Vanrumste, Roberta Grech, Joseph Muscat, Wim De Clercq, Anneleen Vergult, Yves D’Asseler, Kenneth P Camilleri, Simon G Fabri, Sabine Van Huffel, and Ignace Lemahieu.  Review on solving the forward problem in EEG source analysis.  _Journal of NeuroEngineering and Rehabilitation_ , 4(1):46, December 2007.  ISSN 1743-0003.  doi: 10.1186/1743-0003-4-46.  URL https://jneuroengrehab.biomedcentral.com/articles/10.1186/1743-0003-4-46.
  * Halton [1960] J. H. Halton.  On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals.  _Numerische Mathematik_ , 2(1):84–90, December 1960.  ISSN 0029-599X, 0945-3245.  doi: 10.1007/BF01386213.  URL http://link.springer.com/10.1007/BF01386213.
  * Haueisen et al. [1997] J. Haueisen, C. Ramon, M. Eiselt, H. Brauer, and H. Nowak.  Influence of tissue resistivities on neuromagnetic fields and electric potentials studied with a finite element model of the head.  _IEEE Transactions on Biomedical Engineering_ , 44(8):727–735, August 1997.  ISSN 00189294.  doi: 10.1109/10.605429.  URL http://ieeexplore.ieee.org/document/605429/.
  * Haueisen et al. [2002] J. Haueisen, D.S. Tuch, C. Ramon, P.H. Schimpf, V.J. Wedeen, J.S. George, and J.W. Belliveau.  The Influence of Brain Tissue Anisotropy on Human EEG and MEG.  _NeuroImage_ , 15(1):159–166, January 2002.  ISSN 10538119.  doi: 10.1006/nimg.2001.0962.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811901909620.
  * Haueisen et al. [1995] Jens Haueisen, Ceon Ramon, Piotr Czapski, and Michael Eiselt.  On the influence of volume currents and extended sources on neuromagnetic fields: A simulation study.  _Annals of Biomedical Engineering_ , 23(6):728–739, November 1995.  ISSN 0090-6964, 1573-9686.  doi: 10.1007/BF02584472.  URL http://link.springer.com/10.1007/BF02584472.
  * Herman and Usher [2017] Jon Herman and Will Usher.  SALib: An open-source Python library for Sensitivity Analysis.  _Journal of Open Source Software_ , 2(9):97, January 2017.  ISSN 2475-9066.  doi: 10.21105/joss.00097.  URL https://joss.theoj.org/papers/10.21105/joss.00097.
  * Huang et al. [2019] Yu Huang, Abhishek Datta, Marom Bikson, and Lucas C. Parra.  Realistic volumetric-approach to simulate transcranial electric stimulation—ROAST—a fully automated open-source pipeline.  _Journal of Neural Engineering_ , 16(5):056006, July 2019.  ISSN 1741-2552.  doi: 10.1088/1741-2552/ab208d.  URL https://doi.org/10.1088/1741-2552/ab208d.  Publisher: IOP Publishing.
  * Iacono et al. [2015] Maria Ida Iacono, Esra Neufeld, Esther Akinnagbe, Kelsey Bower, Johanna Wolf, Ioannis Vogiatzis Oikonomidis, Deepika Sharma, Bryn Lloyd, Bertram J. Wilm, Michael Wyss, Klaas P. Pruessmann, Andras Jakab, Nikos Makris, Ethan D. Cohen, Niels Kuster, Wolfgang Kainz, and Leonardo M. Angelone.  MIDA: A Multimodal Imaging-Based Detailed Anatomical Model of the Human Head and Neck.  _PLOS ONE_ , 10(4):e0124126, April 2015.  ISSN 1932-6203.  doi: 10.1371/journal.pone.0124126.  URL https://dx.plos.org/10.1371/journal.pone.0124126.
  * Jochmann et al. [2011] Thomas Jochmann, Daniel Güllmar, Jens Haueisen, and Jürgen R. Reichenbach.  Influence of tissue conductivity changes on the EEG signal in the human brain – A simulation study.  _Zeitschrift für Medizinische Physik_ , 21(2):102–112, May 2011.  ISSN 09393889.  doi: 10.1016/j.zemedi.2010.07.004.  URL https://linkinghub.elsevier.com/retrieve/pii/S0939388910000966.
  * Katoch et al. [2019] Nitish Katoch, Bup Kyung Choi, Saurav Z. K. Sajib, EunAh Lee, Hyung Joong Kim, Oh In Kwon, and Eung Je Woo.  Conductivity Tensor Imaging of In Vivo Human Brain and Experimental Validation Using Giant Vesicle Suspension.  _IEEE Transactions on Medical Imaging_ , 38(7):1569–1577, July 2019.  ISSN 1558-254X.  doi: 10.1109/TMI.2018.2884440.  Conference Name: IEEE Transactions on Medical Imaging.
  * Lanfer et al. [2012] B. Lanfer, M. Scherg, M. Dannhauer, T.R. Knösche, M. Burger, and C.H. Wolters.  Influences of skull segmentation inaccuracies on EEG source analysis.  _NeuroImage_ , 62(1):418–431, August 2012.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2012.05.006.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811912004946.
  * Latikka et al. [2001] J.A. Latikka, J.A. Hyttinen, T.A. Kuurne, H.J. Eskola, and J.A. Malmivuo.  The conductivity of brain tissues: comparison of results in vivo and in vitro measurements.  In _2001 Conference Proceedings of the 23rd Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , volume 1, pages 910–912, Istanbul, Turkey, 2001. IEEE.  ISBN 978-0-7803-7211-5.  doi: 10.1109/IEMBS.2001.1019092.  URL http://ieeexplore.ieee.org/document/1019092/.
  * Malmivuo and Plonsey [1995] Jaakko Malmivuo and Robert Plonsey.  _Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields_.  Oxford University Press, October 1995.  ISBN 978-0-19-505823-9.  doi: 10.1093/acprof:oso/9780195058239.001.0001.  URL https://oxford.universitypressscholarship.com/view/10.1093/acprof:oso/9780195058239.001.0001/acprof-9780195058239.
  * McCann et al. [2019] Hannah McCann, Giampaolo Pisano, and Leandro Beltrachini.  Variation in Reported Human Head Tissue Electrical Conductivity Values.  _Brain Topography_ , 32(5):825–858, September 2019.  ISSN 0896-0267, 1573-6792.  doi: 10.1007/s10548-019-00710-2.  URL http://link.springer.com/10.1007/s10548-019-00710-2.
  * Michel and Brunet [2019] Christoph M. Michel and Denis Brunet.  EEG Source Imaging: A Practical Review of the Analysis Steps.  _Frontiers in Neurology_ , 10:325, 2019.  ISSN 1664-2295.  doi: 10.3389/fneur.2019.00325.  URL https://www.frontiersin.org/article/10.3389/fneur.2019.00325.
  * Montes-Restrepo et al. [2014] Victoria Montes-Restrepo, Pieter van Mierlo, Gregor Strobbe, Steven Staelens, Stefaan Vandenberghe, and Hans Hallez.  Influence of Skull Modeling Approaches on EEG Source Localization.  _Brain Topography_ , 27(1):95–111, January 2014\.  ISSN 0896-0267, 1573-6792.  doi: 10.1007/s10548-013-0313-y.  URL http://link.springer.com/10.1007/s10548-013-0313-y.
  * Nielsen et al. [2018] Jesper D. Nielsen, Kristoffer H. Madsen, Oula Puonti, Hartwig R. Siebner, Christian Bauer, Camilla Gøbel Madsen, Guilherme B. Saturnino, and Axel Thielscher.  Automatic skull segmentation from MR images for realistic volume conductor models of the head: Assessment of the state-of-the-art.  _NeuroImage_ , 174:587–598, July 2018.  ISSN 1095-9572.  doi: 10.1016/j.neuroimage.2018.03.001.
  * Nuwer [2018] Marc R. Nuwer.  10-10 electrode system for EEG recording.  _Clinical Neurophysiology_ , 129(5):1103, May 2018\.  ISSN 1388-2457.  doi: 10.1016/j.clinph.2018.01.065.  URL http://www.sciencedirect.com/science/article/pii/S1388245718300907.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Plonsey and Heppner [1967] R. Plonsey and D. B. Heppner.  Considerations of quasi-stationarity in electrophysiological systems.  _The Bulletin of Mathematical Biophysics_ , 29(4):657–664, December 1967.  ISSN 0007-4985.  doi: 10.1007/BF02476917.
  * Puonti et al. [2020] Oula Puonti, Koen Van Leemput, Guilherme B. Saturnino, Hartwig R. Siebner, Kristoffer H. Madsen, and Axel Thielscher.  Accurate and robust whole-head segmentation from magnetic resonance images for individualized head modeling.  _NeuroImage_ , 219:117044, October 2020.  ISSN 1053-8119.  doi: 10.1016/j.neuroimage.2020.117044.  URL https://www.sciencedirect.com/science/article/pii/S1053811920305309.
  * Rasmussen and Williams [2006] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian processes for machine learning_.  Adaptive computation and machine learning. MIT Press, Cambridge, Mass, 2006.  ISBN 978-0-262-18253-9.  OCLC: ocm61285753.
  * Sadleir and Argibay [2007] R. J. Sadleir and A. Argibay.  Modeling Skull Electrical Properties.  _Annals of Biomedical Engineering_ , 35(10):1699–1712, September 2007.  ISSN 0090-6964, 1573-9686.  doi: 10.1007/s10439-007-9343-5.  URL http://link.springer.com/10.1007/s10439-007-9343-5.
  * Sajib et al. [2016] Saurav Z. K. Sajib, Woo Chul Jeong, Eun Jung Kyung, Hyun Bum Kim, Tong In Oh, Hyung Joong Kim, Oh In Kwon, and Eung Je Woo.  Experimental evaluation of electrical conductivity imaging of anisotropic brain tissues using a combination of diffusion tensor imaging and magnetic resonance electrical impedance tomography.  _AIP Advances_ , 6(6):065109, June 2016.  doi: 10.1063/1.4953893.  URL https://aip.scitation.org/doi/full/10.1063/1.4953893.  Publisher: American Institute of Physics.
  * Saltelli [2008] A. Saltelli, editor.  _Global sensitivity analysis: the primer_.  John Wiley, Chichester, England ; Hoboken, NJ, 2008.  ISBN 978-0-470-05997-5.  OCLC: ocn180852094.
  * Saltelli et al. [2010] Andrea Saltelli, Paola Annoni, Ivano Azzini, Francesca Campolongo, Marco Ratto, and Stefano Tarantola.  Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index.  _Computer Physics Communications_ , 181(2):259–270, February 2010.  ISSN 00104655.  doi: 10.1016/j.cpc.2009.09.018.  URL https://linkinghub.elsevier.com/retrieve/pii/S0010465509003087.
  * Saturnino et al. [2019] Guilherme B. Saturnino, Axel Thielscher, Kristoffer H. Madsen, Thomas R. Knösche, and Konstantin Weise.  A principled approach to conductivity uncertainty analysis in electric field calculations.  _NeuroImage_ , 188:821–834, March 2019.  ISSN 1053-8119.  doi: 10.1016/j.neuroimage.2018.12.053.  URL http://www.sciencedirect.com/science/article/pii/S1053811918322031.
  * Schimpf [2007] P. H. Schimpf.  Application of Quasi-Static Magnetic Reciprocity to Finite Element Models of the MEG Lead-Field.  _IEEE Transactions on Biomedical Engineering_ , 54(11):2082–2088, November 2007.  ISSN 1558-2531.  doi: 10.1109/TBME.2007.895112.  Conference Name: IEEE Transactions on Biomedical Engineering.
  * Schirru et al. [2011] Andrea Schirru, Simone Pampuri, Giuseppe De Nicolao, and Sean McLoone.  Efficient Marginal Likelihood Computation for Gaussian Process Regression.  _arXiv:1110.6546 [stat]_ , October 2011.  URL http://arxiv.org/abs/1110.6546.  arXiv: 1110.6546.
  * Smith et al. [2004] Stephen M. Smith, Mark Jenkinson, Mark W. Woolrich, Christian F. Beckmann, Timothy E. J. Behrens, Heidi Johansen-Berg, Peter R. Bannister, Marilena De Luca, Ivana Drobnjak, David E. Flitney, Rami K. Niazy, James Saunders, John Vickers, Yongyue Zhang, Nicola De Stefano, J. Michael Brady, and Paul M. Matthews.  Advances in functional and structural MR image analysis and implementation as FSL.  _NeuroImage_ , 23 Suppl 1:S208–219, 2004.  ISSN 1053-8119.  doi: 10.1016/j.neuroimage.2004.07.051.
  * Sobol [1967] I.M Sobol.  On the distribution of points in a cube and the approximate evaluation of integrals.  _USSR Computational Mathematics and Mathematical Physics_ , 7(4):86–112, January 1967.  ISSN 00415553.  doi: 10.1016/0041-5553(67)90144-9.  URL https://linkinghub.elsevier.com/retrieve/pii/0041555367901449.
  * Sobol [1976] I.M. Sobol.  Uniformly distributed sequences with an additional uniform property.  _USSR Computational Mathematics and Mathematical Physics_ , 16(5):236–242, January 1976.  ISSN 00415553.  doi: 10.1016/0041-5553(76)90154-3.  URL https://linkinghub.elsevier.com/retrieve/pii/0041555376901543.
  * Sobol [2001] I.M Sobol.  Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates.  _Mathematics and Computers in Simulation_ , 55(1-3):271–280, February 2001.  ISSN 03784754.  doi: 10.1016/S0378-4754(00)00270-6.  URL https://linkinghub.elsevier.com/retrieve/pii/S0378475400002706.
  * Taberna et al. [2021] Gaia Amaranta Taberna, Jessica Samogin, and Dante Mantini.  Automated Head Tissue Modelling Based on Structural Magnetic Resonance Images for Electroencephalographic Source Reconstruction.  _Neuroinformatics_ , January 2021.  ISSN 1539-2791, 1559-0089.  doi: 10.1007/s12021-020-09504-5.  URL http://link.springer.com/10.1007/s12021-020-09504-5.
  * Tadel et al. [2011] François Tadel, Sylvain Baillet, John C. Mosher, Dimitrios Pantazis, and Richard M. Leahy.  Brainstorm: A User-Friendly Application for MEG/EEG Analysis.  _Computational Intelligence and Neuroscience_ , 2011:e879716, April 2011.  ISSN 1687-5265.  doi: 10.1155/2011/879716.  URL https://www.hindawi.com/journals/cin/2011/879716/.  Publisher: Hindawi.
  * The CGAL Project [2020] The CGAL Project.  _CGAL User and Reference Manual_.  CGAL Editorial Board, 5.0.3 edition, 2020.  URL https://doc.cgal.org/5.0.3/Manual/packages.html.
  * Thielscher et al. [2015] Axel Thielscher, Andre Antunes, and Guilherme B. Saturnino.  Field modeling for transcranial magnetic stimulation: A useful tool to understand the physiological effects of TMS?  In _2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , pages 222–225, Milan, August 2015. IEEE.  ISBN 978-1-4244-9271-8.  doi: 10.1109/EMBC.2015.7318340.  URL http://ieeexplore.ieee.org/document/7318340/.
  * Tuch et al. [2001] D. S. Tuch, V. J. Wedeen, A. M. Dale, J. S. George, and J. W. Belliveau.  Conductivity tensor mapping of the human brain using diffusion tensor MRI.  _Proceedings of the National Academy of Sciences_ , 98(20):11697–11701, September 2001.  ISSN 0027-8424, 1091-6490.  doi: 10.1073/pnas.171473898.  URL http://www.pnas.org/cgi/doi/10.1073/pnas.171473898.
  * Vallaghe and Clerc [2009] Sylvain Vallaghe and Maureen Clerc.  A Global Sensitivity Analysis of Three- and Four-Layer EEG Conductivity Models.  _IEEE Transactions on Biomedical Engineering_ , 56(4):988–995, April 2009.  ISSN 1558-2531.  doi: 10.1109/TBME.2008.2009315.  Conference Name: IEEE Transactions on Biomedical Engineering.
  * Vorwerk et al. [2018] Johannes Vorwerk, Robert Oostenveld, Maria Carla Piastra, Lilla Magyari, and Carsten H. Wolters.  The FieldTrip-SimBio pipeline for EEG forward solutions.  _BioMedical Engineering OnLine_ , 17(1):37, March 2018.  ISSN 1475-925X.  doi: 10.1186/s12938-018-0463-y.  URL https://doi.org/10.1186/s12938-018-0463-y.
  * Vorwerk et al. [2019] Johannes Vorwerk, Umit Aydin, Carsten Wolters, and Christopher Butson.  Influence of Head Tissue Conductivity Uncertainties on EEG Dipole Reconstruction.  _Frontiers in Neuroscience_ , 13, June 2019.  doi: 10.3389/fnins.2019.00531.
  * Weinstein et al. [2000] David Weinstein, Leonid Zhukov, and Chris Johnson.  Lead Field Basis for FEM Source Localization.  _Annals of Biomedical Engineering_ , 28(9):1059–1065, 2000.  doi: 10.1114/1.1310220.
  * Windhoff et al. [2013] Mirko Windhoff, Alexander Opitz, and Axel Thielscher.  Electric field calculations in brain stimulation based on finite elements: An optimized processing pipeline for the generation and usage of accurate individual head models.  _Human Brain Mapping_ , 34(4):923–935, 2013.  ISSN 1097-0193.  doi: https://doi.org/10.1002/hbm.21479.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/hbm.21479.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/hbm.21479.
  * Wolters et al. [2006] C.H. Wolters, A. Anwander, X. Tricoche, D. Weinstein, M.A. Koch, and R.S. MacLeod.  Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: A simulation and visualization study using high-resolution finite element modeling.  _NeuroImage_ , 30(3):813–826, April 2006.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2005.10.014.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811905007871.
  * Wu et al. [2018] Zhanxiong Wu, Yang Liu, Ming Hong, and Xiaohui Yu.  A review of anisotropic conductivity models of brain white matter based on diffusion tensor imaging.  _Medical & Biological Engineering & Computing_, 56(8):1325–1332, August 2018.  ISSN 0140-0118, 1741-0444.  doi: 10.1007/s11517-018-1845-9.  URL http://link.springer.com/10.1007/s11517-018-1845-9.
  * Ziegler et al. [2014] Erik Ziegler, Sarah L. Chellappa, Giulia Gaggioni, Julien Q.M. Ly, Gilles Vandewalle, Elodie André, Christophe Geuzaine, and Christophe Phillips.  A finite-element reciprocity solution for EEG forward modeling with realistic individual head models.  _NeuroImage_ , 103:542–551, December 2014.  ISSN 10538119.  doi: 10.1016/j.neuroimage.2014.08.056.  URL https://linkinghub.elsevier.com/retrieve/pii/S1053811914007307.