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

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$ mm³ 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.

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].

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.

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.

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)

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)

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.

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.

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 mm³ 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.

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.