Improving EEG Forward Modeling Using High-Resolution Five-Layer BEM-FMM Head Models: Effect on Source Reconstruction Accuracy

Abstract

Electroencephalographic (EEG) source localization is a fundamental tool for clinical diagnoses and brain-computer interfaces. We investigate the impact of model complexity on reconstruction accuracy by comparing the widely used three-layer boundary element method (BEM) as an inverse method against a five-layer BEM accelerated by the fast multipole method (BEM-FMM) and coupled with adaptive mesh refinement (AMR) as forward solver. Modern BEM-FMM with AMR can solve high-resolution multi-tissue models efficiently and accurately. We generated noiseless 256-channel EEG data from 15 subjects in the Connectome Young Adult dataset, using four anatomically relevant dipole positions, three conductivity sets, and two head segmentations; we mapped localization errors across the entire grey matter from 4000 dipole positions. The average location error among our four selected dipoles is ∼5$\mathrm{m}$$\mathrm{m}$ (±2$\mathrm{m}$$\mathrm{m}$) with an orientation error of ∼${12}^{\circ }$ (±$7^{\circ }$). The average source localization error across the entire grey matter is ∼9$\mathrm{m}$$\mathrm{m}$ (±4$\mathrm{m}$$\mathrm{m}$), with a tendency for smaller errors on the occipital lobe. Our findings indicate that while three-layer models are robust under noiseless conditions, substantial localization errors (10–20$\mathrm{m}$$\mathrm{m}$) are common. Therefore, models of five or more layers may be needed for accurate source reconstruction in critical applications involving noisy EEG data.

1. Introduction

EEG source reconstruction, or source localization, involves estimating the source location of neural activity within the brain from EEG-recorded measurements [1]. Applications include pre-surgical screening for the treatment of drug-resistant epilepsy [2], brain-computer interfaces (BCIs) [3], and the analysis of event-related potentials in a wide range of scenarios [4,5,6]. Due to the large susceptibility of EEG to noise and the ill-posedness of the EEG inverse problem [7], it is important to reduce modeling errors in the forward step of source localization pipelines; this can be carried out by distinguishing the relevant regions of different conductivity present in the head [8,9], utilizing accurate values for such conductivities [10], and matching the geometry and features of the head closely using high-resolution models [11]. For example, the geometry of the intracortical compartments can be very important for source reconstruction using source models where the dipole orientation is perpendicular to the cortex [12].

Many excellent open-source software packages implement EEG source localization, including Brainstorm [13], FieldTrip [14], MNE [15], and EEGLab [16]. These BEM implementations typically use three layers extracted from the subject’s MRI image: scalp, outer skull, and inner skull. For most practical applications, the resolution of these layers is limited to less than 10,000 triangles per layer. One of the main reasons for this limitation is that these packages use the classical potential-based formulation of the BEM [17,18,19], which involves dense system matrices and, as such, is unable to compute leadfields using large, high-resolution models. Complicating this issue is the fact that the intracortical compartments typically involve very narrow gaps [20] between the grey and white matter layers, where the sources are placed. This can cause large numerical inaccuracies in the forward computations, which need to be resolved with refinement techniques such as adaptive mesh refinement (AMR), also known as h-refinement [21,22], or the less common p-refinement [23], which consists of adaptively increasing the polynomial order of the local approximations of the variable of interest (potential or charge) on the mesh triangles. Both refinement techniques cause a computational overhead, which may be avoided in three-layer models as long as the sources are far away from the discretized shells.

Modern automated human-head segmentation tools (FreeSurfer, [24]; SPM12, [25]) can provide high-resolution skin, skull, cerebrospinal fluid (CSF), grey matter (GM), and white matter (WM) layers, among other tissues. The charge-based formulation of the BEM, coupled with fast multipole method (FMM) acceleration, or BEM-FMM [19,26], can rapidly compute forward solutions for high-resolution models involving tens of millions of facets, thereby removing the ∼10,000 triangle limitation of conventional potential-based BEM. Therefore, BEM-FMM could be used as the basis of high-resolution multi-layer EEG source reconstruction pipelines, offering a powerful alternative to high-resolution FEM-based forward solvers [27].

To assess the potential improvements in source reconstruction that may arise from the use of BEM-FMM, we consider the following question: How accurate would the three-layer BEM source reconstruction results be if we assume the high-resolution five-layer BEM-FMM models are ground truth? We investigated this in a modeling experiment for the canonical EEG inverse problem of single-dipole fitting.

There are many studies that consider the influence of different uncertainty factors on EEG source localization: head model complexity [9], conductivity uncertainty [10], and white matter anisotropy [28,29], among others. In all these studies, the finite element method (FEM) was used for forward modeling, i.e., for generating electrode voltages given a dipole source in the brain. Here, we use a high-resolution BEM-FMM to compute the forward solution. This is particularly motivated by the following reasons:

• (i). Studies of dipole localization accuracy using BEM [30,31,32] are relatively less common than those using FEM [9,10,28,29,33]. The use of BEM-FMM allows us to carry out more tests and with higher resolution than conventional BEM.

• (ii). BEM-FMM models have been found to be comparable or superior to FEM in terms of speed and accuracy in certain applications; for example, it has been found [34] that zero-order BEM-FMM is comparable in accuracy to second-order FEM for TMS modeling in the concentric spheres model. It was also found in [34] that zero-order BEM-FMM was the highest-accuracy method that they could implement within their computational constraints using a realistic, high-resolution brain model. This leads them to use zero-order BEM-FMM as their ground truth solution.

• (iii). Prescribing arbitrary finite-length dipoles, or point dipoles—which are the most common mathematical model for the simultaneous firing of a large number of neurons [35]—is an easy task in BEM, as the method inherently allows for arbitrary incident fields anywhere in space. On the other hand, modeling dipoles in FEM is much more challenging: several approaches are available [36], but using them would involve an additional error estimation; see also the St. Venant approach used in [10,37].

With regards to (ii), we must mention that one major advantage of FEM over BEM is that FEM can easily model tissue anisotropy. Although it is possible to model anisotropy in BEM [38], most BEM/BEM-FMM implementations do not have this capability at the time of writing.

• Section 2 describes the materials used throughout the paper and the methods used for the analysis of four dipole locations chosen consistently across subjects;

• Section 3 describes a complementary analysis to the one in Section 2, where we generated error maps for source localization across the entire grey matter surface of each subject;

• Section 4 summarizes our results;

• Section 5 includes a brief discussion and interpretation.

In Supplement 1 (Supplementary Tables S1–S14), we include additional tables of results, and in Supplement 2 (Supplementary Figures S1–S17), we include additional subject error maps and figures.

2. Localization Error in Selected Locations Across Varying Subjects and Parameters

Outline of methodology: We set up the following numerical experiment, which is carried out for a total of 15 subjects, three conductivity sets, two mesh segmentations, and four dipole locations:

• (i). We placed a dipole at a location, ${\mathbf{p}}_0$, on the midsurface between the CSF-GM and GM-WM tissue interfaces of our subject with an orientation, ${\mathbf{q}}_0$, normal to the CSF-GM interface.

• (ii). We simulated a single time sample of noiseless EEG data using the charge-based formulation of the boundary element method with fast multipole method acceleration (BEM-FMM) and adaptive mesh refinement (AMR), over a high-resolution five-layer (seven-compartment) head volume conduction model.

• (iii). We used the simulated EEG data to perform source reconstruction with a low-resolution three-layer head model, which is similar to the ones widely used in EEG source reconstruction. By doing so, we found the best fit for location ${\mathbf{p}}_1$, and orientation ${\mathbf{q}}_1$ was provided by the FieldTrip Toolbox’s [14] source localization procedures.

• (iv). We computed the distance $\parallel {\mathbf{p}}_0-{\mathbf{p}}_{\mathbf{1}}{\parallel }_2$ between both locations and the angle between ${\mathbf{q}}_0$ and ${\mathbf{q}}_1$ to measure the error of the fit.

Five-layer models. We used the SimNIBS package [39] to obtain two segmentations from the T1- and T2-weighted MRI data from 15 subjects of the Connectome Young Adult dataset [40]. One segmentation was obtained using the option mri2mesh, which performs segmentation based on FreeSurfer [24], and the other used headreco, which is based on SPM/CAT [25]. The FreeSurfer and Headreco segmentations generate five main layers: skin, skull, CSF, grey matter, and white matter.

Each skin, skull, and CSF mesh from the FreeSurfer segmentations contains approximately 120,000 triangles, whereas the grey and white matter meshes contain circa 250,000 triangles. Headreco segmentations comprise approximately 150,000 triangles for skin and skull, 100,000 triangles for CSF, 300,000 triangles for grey matter, and 350,000 triangles for white matter. So, in total, our high-resolution FreeSurfer models have over 1 M degrees of freedom. This number is about twice as large in the Headreco models.

Both segmentations include two additional tissue meshes; in the case of FreeSurfer, these are the cerebellum and ventricles, and for Headreco, these are the eyes and ventricles. One may argue that this should limit comparability between segmentations to some degree; however, the location of the two additional compartments is distant from most electrode positions. More details on the MRI data acquisition, spatial coregistration, segmentation, and mesh generation for this dataset can be found in [40,41].

We tested three tissue conductivity sets for each subject, which we labeled IT’IS7 [42,43], VWB7 [10,44], and SimNIBS7 [39]. The conductivity values in IT’IS7 and SimNIBS7 were determined experimentally, and the values in VWB7 were estimated from experimental conductivity values [45,46,47,48,49]. We made no attempt to estimate the individual tissue conductivities of each subject, although we mention that this may be possible using the methods of [50]. All conductivity sets for five-layer models can be found in Table 1, and for additional information on the conductivity of living tissues, see [1].

Each of the chosen conductivity sets is optimized for different scenarios: The IT’IS7 set provides conductivity values for frequencies up to 1$\mathrm{M}$$\mathrm{Hz}$, SimNIBS7 specifies conductivities adequate for TMS and TES (see also [51,52]), and VWB7 has been derived from studies on the uncertainty of EEG source localization with respect to conductivity values.

The only set that provides conductivity values for all the layers is IT’IS7 [43]. For VWB7, no standard cerebellum conductivities are given. We chose a 40–60 % weighted average between the conductivities of grey matter and white matter as a conductivity value for the cerebellum. IT’IS7 assigns the same conductivity value as CSF for the eye’s vitreous humor. Hence, we assign the corresponding CSF conductivity for the eye layer in each Headreco conductivity set. We must remark, however, that the conductivity values of these additional layers are mostly irrelevant in the EEG problem since these are located far away from the scalp electrodes.

Three-layer models. From each FreeSurfer seven-compartment head model, we created decimated (downsampled) three-layer models using the software MeshLab [53]. Namely, we created three new meshes labeled SKIN, SKULL, and BRAIN from the high-resolution skin, skull, and CSF FreeSurfer layers of each subject.

To obtain these new meshes, we first apply a screened Poisson surface reconstruction filter [54], followed by a quadratic edge collapse decimation filter with 14,000 triangles as a target. After this is carried out, we apply a Taubin smoothing filter [55].

We used these three-layer meshes to carry out dipole source localization using the FieldTrip toolbox [14]. In other words, we used the three-layer model as our inverse volume conduction model. For every choice of conductivity for the forward five-layer model, we have the corresponding three-layer conductivities labeled IT’IS3, VWB3, and SimNIBS3; see Table 2. Each conductivity set was chosen to be as consistent as possible with its seven-compartment analog.

Source dipole placement. The following are the four selected locations that we sampled consistently across all subjects:

• (i). dip1 — Posterior wall of the central sulcus, somatosensory cortex, and tangential dipole; Figure 1;

• (ii). dip2 — ${\mathrm{M}1}_{\mathrm{HAND}}$ region, primary motor cortex, and radial dipole; Figure 2;

• (iii). dip3 — Temporal lobe, along the Heschl’s gyri or transverse temporal gyri; Figure 3;

• (iv). dip4 — Medio-temporal region; Figure 4.

The position dip3 is at the superior part of the temporal cortex and has been chosen to be studied in both auditory and language areas located at the posterior or anterior part around the same depth. The position dip4 is relevant in epileptogenic analyses [56,57].

Electrode Placement: We placed 256 electrodes in each subject’s FreeSurfer, Headreco, and decimated FreeSurfer skin meshes; see Figure 5. The high-density 256 dry electrode system that we used is described in [58]; see [59]. Since our dataset [41] does not include data on the three anatomical landmarks of the nasion, left preauricular point (LPA), or right preauricular point (RPA), we did a semi-manual electrode placement on each subject’s skin layer using the FieldTrip function ft_electrode_placement, first using the interactive option, followed by the projection option. In this way, each electrode was assigned to a triangle in the skin layer. We could then retrieve the electrode voltages from the BEM-FMM calculation of the skin potential. This was carried out by taking the value of the potential at the mesh triangle of each electrode.

Forward solution with BEM-FMM: We used the charge-based formulation of the BEM [60,61], together with FMM acceleration [62].

In this formulation of the BEM forward problem, the total electric field, $\mathbf{E}$, is separated into two components, an impressed electric field, ${\mathbf{E}}^i$, and a secondary electric field ${\mathbf{E}}^s$, so that $\mathbf{E}={\mathbf{E}}^i+{\mathbf{E}}^s$. In our setting, the impressed field is created by a current dipole [63,64] in an infinite homogeneous medium having the conductivity of grey matter. Such a current dipole is given by an isotropic current source of strength $+I_0$ and an isotropic current sink of strength $-I_0$, both measured in amperes. For point current dipoles, the length of the displacement vector, $\mathbf{d}$, from the sink to the source is small with respect to the distance of the dipole center, ${\mathbf{r}}_0$, to the observation point, $\mathbf{r}$. In this case, the impressed electric field can be approximated as ${\mathbf{E}}^i=-\nabla {\phi }^i$, where

(1)${\phi }^i\left(\mathbf{r}\right)={\displaystyle \frac{I_0}{4\pi {\sigma }_{\mathtt{GM}}}}{\displaystyle \frac{\mathbf{d}·(\mathbf{r}-{\mathbf{r}}_0)}{|\mathbf{r}-{\mathbf{r}}_0{|}^3}}.$

Here, the constant ${\sigma }_{\mathtt{GM}}$ is the macroscopic conductivity of grey matter tissue. When this constant is measured in ${\mathrm{S} \mathrm{m}}^{-1}={\mathrm{A} \mathrm{V}}^{-1}{\mathrm{m}}^{-1}$}, and the distances are measured in meters, then the potential ${\phi }^i$ is measured in volts as desired.

The secondary electric field is always conservative and arises from the electric charges deposited on each surface of discontinuity, S, of the conductivity upon activation of the impressed electric field. To compute the secondary field, ${\mathbf{E}}^s$, we need to calculate the surface charge density, $\rho$, on S, measured in $\mathrm{C} {\mathrm{m}}^{-2}$ An integral equation can be derived [19] using the quasi-static approximation of Maxwell’s equations and the following boundary conditions:

(2)$\left\{\begin{array}{cc}{\sigma }_{+}{\mathbf{E}}_{+}·\mathbf{n}={\sigma }_{-}{\mathbf{E}}_{-}·\mathbf{n},\hfill & \mathrm{on}S\hfill \\ {\mathbf{E}}_{-}^i={\mathbf{E}}_{+}^i,\hfill & \mathrm{on}S\hfill \end{array}\right..$

Here, $\mathbf{n}\left(\mathbf{r}\right)$ denotes the outward unit normal vector to the point $\mathbf{r}\in S$, and $f_{\pm }\left(\mathbf{r}\right)={lim}_{ϵ\to 0}f(\mathbf{r}\pm ϵ\mathbf{n}\left(\mathbf{r}\right))$ for all $\mathbf{r}\in S$. The function $\sigma \left(\mathbf{r}\right)$ is the conductivity of our head model at an arbitrary point $\mathbf{r}$ in 3D space, measured in ${\mathrm{S} \mathrm{m}}^{-1}$. The resulting integral equation for $\rho$ is given below [26,65,66]:

(3)$\begin{array}{c}\hfill {\displaystyle \frac{\rho \left(\mathbf{r}\right)}{2{\epsilon }_0}}-K\left(\mathbf{r}\right)\mathbf{n}\left(\mathbf{r}\right)·{\int }_S{\displaystyle \frac{\rho \left({\mathbf{r}}^{\prime }\right)}{4\pi {\epsilon }_0}}{\displaystyle \frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^3}}ds\left({\mathbf{r}}^{\prime }\right)\\ \hfill =K\left(\mathbf{r}\right){\mathbf{E}}^i\left(\mathbf{r}\right)·\mathbf{n}\left(\mathbf{r}\right),\end{array}$

where ${\epsilon }_0\approx 8.8541 \times {10}^{-12} {\mathrm{F} \mathrm{m}}^{-1}$ is the permittivity of free space, $\mathbf{n}\left(\mathbf{r}\right)$ is the outward unit normal vector to S at $\mathbf{r}$, and $K:=({\sigma }_{-}-{\sigma }_{+})/({\sigma }_{-}+{\sigma }_{+})$ is the conductivity contrast at $\mathbf{r}\in S$. Since the surface element $ds\left({\mathbf{r}}^{\prime }\right)$ is measured in ${\mathrm{m}}^2$, and both the unit normal $\mathbf{n}$ and the conductivity contrast K are dimensionless quantities, all terms in Equation (3) are measured in ${\mathrm{V} \mathrm{m}}^{-1}$.

Equation (3) can be discretized using the Galerkin method. For a triangular mesh model comprising M facets, which approximates the geometry of each surface of discontinuity S, the discrete version of Equation (3) at the m-th facet reads

(4)$\begin{array}{c}\hfill {\displaystyle \frac{{\rho }_m}{2{\epsilon }_0}}-K_m{\mathbf{n}}_m·\sum _{n=1;n\ne m}^M{\displaystyle \frac{{\rho }_n}{4\pi {\epsilon }_0}}{\displaystyle \frac{{\mathbf{r}}_m-{\mathbf{r}}_n}{|{\mathbf{r}}_m-{\mathbf{r}}_n{|}^3}}\\ \hfill =K_m{\mathbf{E}}^i\left({\mathbf{r}}_m\right)·{\mathbf{n}}_m.\end{array}$

Here, ${\rho }_m$, $K_m$, ${\mathbf{n}}_m$, and ${\mathbf{r}}_m$ are the charge density, conductivity contrast, outward unit normal, and triangle center at the m-th facet, respectively. The summation term is an n-body computation that we accelerated using FMM [62,67,68]. In our computations, we used the implementation of the fast multipole method provided by the FMM3D library [69]. The FMM-accelerated summation was computed to a digit precision of $1 \times {10}^{-6}$. For reference, the typical FMM digit precision in most applications ranges from $1 \times {10}^{-4}$ to $1 \times {10}^{-2}$.

To estimate ${\rho }_m$ at every facet, successive approximations, ${\rho }_m^{\left(0\right)},{\rho }_m^{\left(1\right)},\cdots ,{\rho }_m^{\left(N\right)}$, of the charge density are computed using FMM acceleration. After each step, we applied the generalized residual method (GMRES) to find an approximation to ${\rho }_m$ of the form

(5)${\rho }_m={\alpha }_0{\rho }_m^{\left(0\right)}+\cdots +{\alpha }_N{\rho }_m^{\left(N\right)},$

where the scalars ${\alpha }_i$ are chosen to minimize the residual:

(6)${\left|{\displaystyle \frac{{\rho }_m}{2{\epsilon }_0}}-K_m{\mathbf{n}}_m·\left(\sum _{\begin{array}{c}n=1\\ n\ne m\end{array}}^M{\displaystyle \frac{{\rho }_n}{4\pi {\epsilon }_0}}{\displaystyle \frac{{\mathbf{r}}_m-{\mathbf{r}}_n}{|{\mathbf{r}}_n-{\mathbf{r}}_m{|}^3}}+{\mathbf{E}}^i\left({\mathbf{r}}_m\right)\right)\right|}^2.$

See [19] for more details on the computation of EEG forward solutions using BEM-FMM. In our case, we stopped iterating when the GMRES residual was below $1 \times {10}^{-5} {\mathrm{V}}^2 {\mathrm{m}}^{-2}$. For reference, typical BEM-FMM GMRES residual thresholds range from $1 \times {10}^{-4}$ to $1 \times {10}^{-3} {\mathrm{V}}^2 {\mathrm{m}}^{-2}$; see [21] for further details on the precision of BEM-FMM utilized.

Adaptive mesh refinement: Adaptive mesh refinement (AMR) is a general technique that ensures the convergence of the discretized solution to the BEM problem to the true analytic solution [70]. The AMR method in the context of BEM-FMM was introduced in [22] (also see [21]).

AMR consists of subdividing mesh triangles according to a cost function until a stopping criterion, based on the convergence of the solution over a region of interest, is reached. In our case, the cost function for the m-th triangular facet is $C_m=\left|{\rho }_m\right|·A_m$, where ${\rho }_m$ is the charge density in the facet, and $A_m$ is its area. In other words, our cost function is the total charge magnitude in each triangle. The triangles with the top 1% cost over all triangles are subdivided into four congruent triangles. After this, the surface charges are recomputed using BEM-FMM, new costs are computed, and a new AMR step begins. The AMR process is terminated whenever the change in electrode voltages is less than 1% relative to the previous AMR step; to be precise, let $V_n\in {\mathbb{R}}^{256}$ be the measured electrode potentials at the n-th AMR step. We stop the iterative procedure at the n-th step provided that

(7)$\parallel V_n-V_{n-1}{\parallel }_2/{\parallel V_{n-1}\parallel }_2<0.01,$

where ${\parallel ·\parallel }_2$ is the 2-norm in ${\mathbb{R}}^{256}$. This stopping criterion has a simple interpretation: Once the relative change in voltage values is very small (less than 1%), we consider that the AMR algorithm has converged to a good approximation of the analytical solution.

When we perform AMR, all layers of the model are refined except for the skin layer, which is our region of interest, and it is used to compare the relative change from one AMR step to the next.

The AMR implementation that we used is the one utilized and described in [21]. In [71], the reader can find error analyses of BEM-FMM with AMR against the analytical solution of the concentric spheres model for EEG. We did not perform particular error estimates against the true solution of the EEG forward model in the five-layer model. We mention, however, that the AMR solution is guaranteed to converge to the analytical solution under certain assumptions [70].

AMR has been incorporated to ensure the accuracy of the five-layer forward solution since the proximity of the dipole source singularity to the triangular elements of the mesh can cause large numerical errors; see Figure 6. AMR has never been applied to the dipole fitting with three-layer models, as this is not part of Fieldtrip’s three-layer bemcp BEM implementation.

EEG source localization: For EEG source localization, we used the FieldTrip MATLAB Toolbox [14]. We used the functions ft_prepare_headmodel with the three-layer decimated meshes, using each of the conductivity values in Table 2. A single-dipole source localization (with free dipole position and orientation) was carried out using nonlinear optimization [72] via the FieldTrip function ft_dipolefitting. The forward solver used in the dipole-fitting procedure was bemcp [73], which implements the classical surface potential formulation of the EEG forward problem [17,74] using a three-layer BEM model. We note that this forward solution does not include AMR.

For each subject and dipole position, the input data for the source localization were the data generated from the five-layer BEM-FMM models, using each of the five-layer conductivity sets: IT’IS7, SimNIBS7, and VWB7. The subsequent dipole localization was carried out using the corresponding three-layer conductivities: IT’IS3, SimNIBS3, and VWB3.

The “best fit” for the dipole was reported according to a least-squares approximation of the ground truth EEG data by the EEG data generated using the fitted dipole. The accuracy of the dipole fit is, thus, measured in terms of residual variance (RV), which is the variance in EEG data unexplained by the dipole fit. Therefore, smaller values of residual variance indicate a better fit.

We carried out two dipole-fitting strategies, and in every case, we reported the predicted dipole that has the smallest residual variance. In the first method, we start a nonlinear fit with the true source dipole position as an initial value for the inverse problem [72]. In the second method, we prepare a grid with a resolution of 5 using the function ft_prepare_sourcemodel, and a grid search is performed from every grid position within the skull. The best initial position in the grid is then used as the initial point of a nonlinear fit.

3. Localization Error Maps on the Grey Matter Surface of Each Subject

Outline of the procedure: For each subject, we placed 4000 point dipoles along the midsurface between grey and white matter shells and computed the corresponding forward solution using BEM-FMM. These 4000 forward solutions were computed with five-layer Headreco segmentation models and conductivity values given by the SimNIBS7 set; Table 1. Due to the large number of dipoles, we used a refinement heuristic called “b-refinement”, which was recently introduced in [71] and dramatically reduces the number of iterations required to obtain an accurate forward EEG solution.

The dipole fitting is carried out using the three-layer BEM model with SimNIBS3 conductivities, as in Section 2, with the only difference being that this time, we do not perform the computationally expensive grid search. We observed that in most cases, the distance between the fitted dipoles in the grid search and in the search starting at the source point is of the order of 10 $\mathrm{\mu }$$\mathrm{m}$ or less. Therefore, the error calculations should not differ considerably between both methods.

Once the 4000 dipole fits have been calculated, we compute the distance between each fitted dipole and the dipole source. To produce an image of the distance error over the entire grey matter, we interpolate these results according to the reciprocal of the distance of each triangle center to each dipole source within a given interpolation radius, R. In other words, to a triangle with center $\mathbf{r}$, we assign the error value

(8)${\displaystyle \frac{{\sum }_{i=1}^k{ϵ}_i/|\mathbf{r}-{\mathbf{s}}_i|}{{\sum }_{i=1}^{k}1/|\mathbf{r}-{\mathbf{s}}_i|}},$

where $\{{\mathbf{s}}_1,\cdots ,{\mathbf{s}}_k\}$ is the set of dipole sources at a distance of $\le R$ from $\mathbf{r}$, and ${ϵ}_i$ is the fitting distance error at the dipole with center ${\mathbf{s}}_i$. In Figure 7 and Figure 8, we exhibit some of the two error maps we generated.

The method of b-refinement: In conventional AMR, the cost function $C_m^{\left(n\right)}$ at step n is given by ${\rho }_m^{(n-1)}$, i.e., the surface charge density computed in the previous AMR step. In b-refinement (or b-AMR), instead of computing the full solution to the charge density in the previous step, we only compute its zero-order approximation, which is given by the incident field at the mesh triangle centers. This can be expressed as

(9)$C_m=\left|{\mathbf{E}}_m^i\right|·A_m,$

where ${\mathbf{E}}^i$ is the incident field at the m-th facet, and $A_m$ is the area of the m-th facet. Equation (9) does not explicitly show its dependence on the refinement step n, but it varies with it since the meshes are subdivided according to the previous computation of $C_m$. In [71], the authors study the accuracy of b-refinement using a concentric spheres model, where the analytic solution is known, as well as a high-resolution five-layer head model using BEM-FMM with conventional AMR as a reference solution. Their findings show that in the high-resolution headmodel, b-refinement exhibits an RDM error of under 2% and a 2-norm error of under 5% for the EEG forward solution.

4. Results

Summary of results: The average distance error in source localization for our four selected locations (across all subjects, segmentations, and conductivities) caused by the discrepancy in charge-based BEM-FMM and potential-based BEM is ∼1$\mathrm{m}$$\mathrm{m}$ in a three-layer model. The difference in the forward methods, together with the discrepancy between the three- and five-layer models, causes an average distance error of ∼4.5$\mathrm{m}$$\mathrm{m}$.

The average error over the entire grey matter surface —using the SimNIBS7 conductivity set, and the Headreco segmentation—is ∼9$\mathrm{m}$$\mathrm{m}$ across all subjects. The patterns of regions of high- and low-dipole-fitting quality appear to be unique to every subject. The main trend that we observed is that the quality of dipole fit tends to be better in the lateral and posterior regions of the brain for most of the subjects, although there are clear exceptions to this.

Three-layer model comparison: FieldTrip does not use the charge-based BEM-FMM in its forward engine; instead, it uses the classical potential-based formulation of the BEM [17,74]. Because of this, we made a model comparison to estimate the baseline accuracy of FieldTrip’s potential-based three-layer inverse model (which does not use AMR) under the assumption that the forward charge-based BEM-FMM+AMR solution calculated using the three-layer model is ground truth. This is a necessary step since there is no analytic solution available, and by knowing the discrepancy due to model differences, we will be able to better assess the decrease in the accuracy of dipole fitting due to the inclusion of additional layers.

First, we calculated the potentials in the entire SKIN layer using a forward three-layer model with charge-based BEM-FMM and AMR. Then, we retrieved readings for the electrode potentials from this calculation using the potential value at the SKIN triangle corresponding to each electrode position. With these electrode voltages, we generated a raw EEG dataset (compatible with FieldTrip software) consisting of a single time sample. Finally, we performed EEG source reconstruction from the simulated EEG using a FieldTrip three-layer volume conduction model with the same conductivity values used in the forward model. See Figure S17 for a diagrammatic representation of the three-layer model comparison process.

The FieldTrip inverse model retrieves the source dipole position from three-shell BEM-FMM+AMR data with an average distance error of ∼1$\mathrm{m}$$\mathrm{m}$ (over all dipoles, subjects, and models) and an average angle error of ${0.75}^{\circ }$. Detailed tables with three-layer model comparison values can be found in Table S13. Comparison tables without the use of AMR for the input data can be found in Table S14.

Five-layer vs. three-layer results: We generated the forward BEM-FMM+AMR solution from each dipole location and each five-layer model, assuming that these solutions are ground truth. Then, we performed dipole fitting using the three-layer model of the corresponding conductivity as an inverse model (e.g., when using IT’IS7 as the conductivity set for our forward computations, we used IT’IS3 as the conductivity set for the inverse computation); see Figure S16 for a flowchart figure. We measured the distance error and orientation error in each case.

We obtained an average distance error (across all dipole positions and subjects) of ∼4$\mathrm{m}$$\mathrm{m}$ (∼1$\mathrm{m}$$\mathrm{m}$ std.) using FreeSurfer meshes and an average distance error of ∼4$\mathrm{m}$$\mathrm{m}$ (∼2.5$\mathrm{m}$$\mathrm{m}$ std.) using Headreco meshes. The error in dipole orientation is more pronounced, with ∼${15}^{\circ }$ average (∼${10}^{\circ }$ std.) for FreeSurfer models and ∼${20}^{\circ }$ average (∼${20}^{\circ }$ std.) for Headreco models. See Table 3 for dipole-wise results.

For additional tables containing data related to the 3-layer model comparison see Tables S5, S6, S13, and S14. For model-wise results for three-layer vs. five-layer models with and without AMR, see Tables S9–S12.

Subject error maps: We located 4000 dipoles in the midsurface between the grey and white matter of every subject. This placement was carried out by first creating a k-means clustering of the grey matter shell and then placing dipoles in the closest midsurface point to each cluster center. We generated the forward solution for each of the 4000 dipoles in isolation using BEM-FMM coupled with b-AMR, and we performed dipole fitting on each of the corresponding voltage readings. This produced a data vector of 4000 distance error values corresponding to each of the dipole centers. Once we had this data, we interpolated it to obtain error values for all the triangle centers of the grey matter mesh; see Figure 7 and Figure 8. The complete set of errormap figures can be found in Figures S1–S15.

Table 4 summarizes the basic statistics associated with each of the error maps that we produced. Subject 122317 is an outlier of our sample; for this subject, the left hemisphere presents much larger localization errors than the right hemisphere. All other patients present largely symmetrical and smooth error maps, and most of them exhibit better accuracy in the posterior and lateral areas of the grey matter.

5. Discussion

Selected dipole fittings: Among the four selected dipoles considered, dip4 has the worst localization results, which is perhaps to be expected since its location in the medio-temporal region (Figure 4) is the deepest of the four. More surprisingly, the position of dip2 (a radial dipole located in the ${\mathrm{M}1}_{\mathrm{HAND}}$ region) was the one that gave the second-to-worst results. The reason may be that the three-layer models do not include CSF, and the closer proximity of the dipole to the highly conductive CSF layer creates a larger discrepancy between the five- and three-layer models. However, one should not over-interpret these results, as more tests for radial dipoles with varying thicknesses of neighboring CSF layers may be needed before a conclusion can be made.

Advantages and shortcomings of the three-layer model: The residual variances we obtained (circa 0.0015 on average) are extremely small, which indicates that the nonlinear dipole fitting using the three-layer model explains most of the EEG data variance with the predicted positions. This indicates that these are good-quality fits. The three-layer model retrieves dipole locations and orientations with good average accuracy from noiseless data. However, when inspecting the histogram of the dipole-fitting errors, one can see common occurrences of errors above 2 cm (Figure 9). Given that we used AMR and b-AMR for our experiments, these errors can be attributed mostly to the discrepancy between model complexity. This suggests that more complex models of five or more layers may be necessary for accurate source reconstruction.

The reader should be aware that, although we did not attempt to model noise in this study, the dipole-fitting errors may be significantly larger when the data is noisy, as is the case with real EEG data. Despite its shortcomings, the three-layer BEM model remains a useful tool for source localization, as it is very computationally efficient and does not typically require refinement methods.

Advantages and challenges of the five-layer model: Inclusion of the grey and white matter tissues is a major factor in the variation of dipole orientation (see Figure 10); this is in agreement with the findings in [10].

Model-wise, the conductivity set that produces the best agreement between the three-layer and five-layer models is IT’IS. A possible explanation is that this set has the lowest conductivity contrast between grey matter and CSF. We also observe (comparing Tables S9 and S10 against Tables S11 and S12) that the dipole fits for the forward data generated with Headreco models have higher variability than those using the data from the FreeSurfer models (both in terms of location and orientation). The use of AMR theoretically eliminates the discrepancies caused by the different number of triangles in the initial layers; therefore, a possible explanation for the different variabilities is that Headreco meshes use a more realistic skull layer than the FreeSurfer meshes. This explanation aligns with the results in [11]. However, a formal assessment of the effect of mesh uncertainty would require a more detailed analysis using techniques similar to those in [10,44].

Figure 6 shows that it is necessary to incorporate some form of mesh refinement to have accurate forward computations and, in turn, an accurate dipole fitting based on five-layer BEM models. Certain dipole positions require up to 12 conventional AMR steps before they can reduce the electrode error sufficiently. For this reason, a practical implementation of dipole fitting using BEM-FMM would require the use of b-AMR or some other accurate heuristic. This is in contrast to three-layer models: the lack of intracortical compartments implies that the dipolar sources will be typically well-separated from the model meshes, and thus, the error caused by the lack of AMR is much lower.

Implications for practical source reconstruction: The error maps in Figure 7 and Figure 8 suggest that single-dipole EEG source localization may be improved for sources located on the frontal and parietal lobes with the use of higher-resolution models, as these are the areas where the three-layer BEM model typically performs worse. However, this is not a general rule, as for certain subjects, the three-layer BEM performs worse on the occipital lobe; see for example Figures S1 and S9.

The method of b-AMR with BEM-FMM appears to be feasible for practical source reconstruction with five-layer models, given our results and the time-performance and error analysis results presented in [71]. Given that our proposed method is also applicable to three- and four-layer models, it would be interesting to analyze the possible performance improvements with high-resolution versions of these mesh models. A future direction to better assess the source localization improvements of high-resolution models would be an analysis using phantoms [12] or experimental subject data [75].

Author Contributions

G.N.P. conceptualization, methodology, software, validation, formal analysis, investigation, visualization, writing—original draft, writing—review & editing; W.A.W. methodology, investigation, software, writing—review & editing; R.C.M. software, validation, formal analysis, investigation, visualization, data curation, writing—review & editing; P.L. software, validation, formal analysis, visualization, data curation, writing—review & editing; J.H. conceptualization, data curation, supervision, writing—review & editing; B.M. conceptualization, supervision, writing—review & editing; T.R.K. conceptualization, supervision, writing—review & editing; K.W. conceptualization, supervision, writing—review & editing; G.M.N. conceptualization, methodology, data curation, supervision, writing—review & editing; T.R. conceptualization, supervision, writing—review & editing; S.N.M. conceptualization, methodology, software, supervision, project administration, funding acquisition, writing—review & editing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The code used to generate the results in this paper is available in the following repository: https://github.com/guillermo-nponasso/EEG-BEMFMM-5Layer.

Acknowledgments

The authors would like to thank Carsten Wolters for his help regarding conductivity values for the three-layer model.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Dipole placement in the transverse, coronal, and sagittal planes, respectively, for dip1 (posterior wall of the central sulcus) in subject 110411. The maroon dot indicates the dipole position on each view. The dashed line indicates the dipole orientation, and the arrow pointing towards the dipole dot is oriented perpendicularly to the dashed line.

Figure 2. Dipole placement in the transverse, coronal, and sagittal planes, respectively, for dip2 ([Forumla omitted. See PDF.]) in subject 130013. The maroon dot indicates the dipole position on each view. The dashed line indicates the dipole orientation, and the arrow pointing towards the dipole dot is oriented perpendicularly to the dashed line.

Figure 3. Dipole placement in the transverse, coronal, and sagittal planes, respectively, for dip3 (temporal lobe and Heschl’s gyrus) in subject 117122. The maroon dot indicates the dipole position on each view. The dashed line indicates the dipole orientation, and the arrow pointing towards the dipole dot is oriented perpendicularly to the dashed line.

Figure 4. Dipole placement in the transverse, coronal, and sagittal planes, respectively, for dip4 (medio-temporal lobe) in subject 122620. The maroon dot indicates the dipole position on each view. The dashed line indicates the dipole orientation, and the arrow pointing towards the dipole dot is oriented perpendicularly to the dashed line.

Figure 5. High-density EEG system of 256 electrodes placed on the skin layer of Connectome subject 110411.

Figure 6. A comparison between the five-layer EEG and MEG forward solutions for dip1 and subject 110411 using the SimNIBS7 conductivity set with and without AMR. In Figure (a), we see the sagittal view of the subject. The colors of the lines correspond to the following layers: white (outer), skin; teal, skull; orange, CSF; apricot, GM; and white (inner) WM. In Figures (b,c), the colors indicate the value of the potential differences in [Forumla omitted. See PDF.] with respect to a zeroth electrode placed on the forehead of the subject. In Figures (d,e), the values represent the strength of the magnetic flux in [Forumla omitted. See PDF.]10[Forumla omitted. See PDF.][Forumla omitted. See PDF.] away from the subject’s skin. Notice the vast differences in the EEG solution profiles: without AMR (b), the forward solution does not exhibit the expected skin voltage profile of a tangential dipole; instead, it shows a profile that would be expected from a radial dipole. However, when AMR is used (c), we observe the typical profile pattern of a tangential dipole. In this particular instance, the source location error with the three-layer BEM for the potentials in (c) is 5.6[Forumla omitted. See PDF.][Forumla omitted. See PDF.] with an orientation error of [Forumla omitted. See PDF.], and the location error for the potentials in (b) is 33[Forumla omitted. See PDF.][Forumla omitted. See PDF.] with an orientation error of [Forumla omitted. See PDF.]. Discrepancies in the profile are also observable in the MEG forward solutions (d,e), although somewhat less pronounced.

Figure 7. Error map for subject 117122: the lighter color corresponds to a better dipole fit location accuracy when using the three-layer BEM with five-layer BEM-FMM+b-AMR as the forward method. (a) transverse view; (b) left sagittal view; (c) anterior coronal view; (d) posterior coronal view.

Figure 8. Error map for subject 120111. (a) transverse view; (b) left sagittal view; (c) anterior coronal view; (d) posterior coronal view.

Figure 9. Histogram of dipole-fitting error ([Forumla omitted. See PDF.] over 4000 well-separated dipole locations across the grey matter surface of subject 117122.

Figure 10. Volumetric plots using a logarithmic scale of the electric potential ([Forumla omitted. See PDF.] around the source dipole dip1 (left hemisphere, posterior wall of the central sulcus opposite to M1 hand region) in subject 110411 for transverse (top anterior; bottom posterior), coronal (top superior; bottom inferior), and sagittal (left posterior; right anterior) views. The colors indicate the varying values of the potential, with the red shades indicating positive values and the blue shades negative values. In the absence of any tissues, the volumetric plots would exhibit mirror symmetry across the plane passing through the dipole center and perpendicular to the dipole direction. The presence of the white matter (highlighted in teal) and the grey matter (highlighted in apricot color) produce a symmetry break that would not occur in a three-layer model; this is particularly noticeable in the sagittal view. Clearly, the symmetry break will contribute to errors in the determination of dipole orientation.

Table 1

Conductivity values (${\mathrm{S} \mathrm{m}}^{-1}$) used in the five-layer/seven-compartment models.

	VWB7	IT’IS7	SimNIBS7
Skin	0.430	0.147	0.465
Skull	0.010	0.0179	0.010
CSF	1.790	1.880	1.654
GM	0.330	0.419	0.275
WM	0.140	0.348	0.126
Cerebellum	0.216	0.577	0.126
Ventricles	1.790	1.880	1.654
Eyes	1.790	1.880	1.654

Table 2

Conductivity values (${\mathrm{S} \mathrm{m}}^{-1}$) for the three-layer models. The values for SKIN and SKULL are taken to be the same as in the corresponding tissues of Table 1. The conductivity of the combined BRAIN layer in the three-layer model has been taken to be 0.33 ${\mathrm{S} \mathrm{m}}^{-1}$ according to [45], except for IT’IS3, where their own value was provided in the dataset.

	VWB3	IT’IS3	SimNIBS3
SKIN	0.430	0.147	0.465
SKULL	0.010	0.0179	0.010
BRAIN	0.330	0.375	0.330

Table 3

Dipole- fitting errors per dipole position over all five-layer models using a BEM-FMM + AMR forward solver for the input voltages and a three-layer classical BEM solver for dipole fitting. The label RV denotes the residual variance of the fit, and AMR denotes the total number of AMR steps required for the convergence of the forward solution.

Dipole	mm	deg	RV	AMR
dip1	4.13	8.71	$6.88 \times {10}^{-4}$	6.22
dip2	5.01	14.62	$8.64 \times {10}^{-4}$	5.66
dip3	2.43	14.04	$9.87 \times {10}^{-4}$	4.59
dip4	6.42	11.00	$8.39 \times {10}^{-4}$	6.23
(a) Averages
dip1	1.65	4.66	$6.39 \times {10}^{-4}$	3.23
dip2	1.85	8.86	$5.37 \times {10}^{-4}$	2.08
dip3	1.13	7.14	$5.69 \times {10}^{-4}$	2.51
dip4	3.87	5.75	$5.89 \times {10}^{-4}$	4.14
(b) Standard Deviations

Table 4

Table of average dipole-fitting errors and standard deviations for 15 subjects of the Connectome Young Adult dataset. Here, dipoles located at the cerebellum and brain stem have been excluded from the computation of averages and standard deviations. Subject 122317 is an outlier in this table.

Subject	Avg. Error (mm)	St. Dev. (mm)
110411	10.28	3.63
117122	9.73	4.11
120111	10.11	4.48
122317	17.71	11.57
122620	5.99	2.27
124422	11.15	2.85
128632	8.05	3.47
130013	7.11	2.44
131722	9.71	3.75
138534	10.17	2.76
149337	7.41	2.51
149539	6.05	2.43
151627	5.57	2.29
160123	7.08	3.03
198451	7.41	3.07

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/bioengineering11111071/s1, Table S1: AMR—FreeSurfer; Table S2: non-AMR—FreeSurfer; Table S3: AMR—Headreco; Table S4: non-AMR—Headreco; Table S5: AMR results per dipole—3-layer model comparison; Table S6: Non-AMR Averages per dipole—3-layer model comparison; Table S7: AMR—all 7-shell models; Table S8: non-AMR results per dipole—all 7-shell models; Table S9: AMR—FreeSurfer results per model and dipole; Table S10: non-AMR—FreeSurfer results per model and dipole; Table S11: AMR—Headreco results per model and dipole; Table S12: non-AMR—Headreco results per model and dipole; Table S13: AMR—3-layer model comparison results per model and dipole; Table S14: non-AMR—3-layer model comparison results per model and dipole. Figure S1: Error map for subject 110411; Figure S2: Error map for subject 117122; Figure S3: Error map for subject 120111; Figure S4: Error map for subject 122317; Figure S5: Error map for subject 122620; Figure S6: Error map for subject 124422; Figure S7: Error map for subject 128632; Figure S8: Error map for subject 130013; Figure S9: Error map for subject 131722; Figure S10: Error map for subject 138534; Figure S11: Error map for subject 149337; Figure S12: Error map for subject 149539; Figure S13: Error map for subject 151627; Figure S14: Error map for subject 160123; Figure S15: Error map for subject 198451; Figure S16: A flowchart for the 3-layer model comparison process; Figure S17: A flowchart for the performance study of 3-layer inverse models using 5-layer forward models.

References

1. Knösche, T.R.; Haueisen, J. EEG/MEG Source Reconstruction; Springer: Cham, Switzerland, 2022.

2. van Mierlo, P.; Vorderwülbecke, B.J.; Staljanssens, W.; Seeck, M.; Vulliémoz, S. Ictal eeg source localization in focal epilepsy: Review and future perspectives. Clin. Neurophysiol.; 2020; 131, pp. 2600-2616. [DOI: https://dx.doi.org/10.1016/j.clinph.2020.08.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32927216]

3. Guttmann-Flury, E.; Sheng, X.; Zhu, X. Channel selection from source localization: A review of four EEG-based brain–computer interfaces paradigms. Behav. Res.; 2023; 55, pp. 1980-2003. [DOI: https://dx.doi.org/10.3758/s13428-022-01897-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35794417]

4. Foti, D.; Weinberg, A.; Dien, J.; Hajcak, G. Event-related potential activity in the basal ganglia differentiates rewards from nonrewards: Temporospatial principal components analysis and source localization of the feedback negativity. Hum. Brain Mapp.; 2011; 32, pp. 2207-2216. [DOI: https://dx.doi.org/10.1002/hbm.21182] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21305664]

5. Hämäläinen, J.A.; Lohvansuu, K.; Ervast, L.; Leppänen, P.H.T. Event-related potentials to tones show differences between children with multiple risk factors for dyslexia and control children before the onset of formal reading instruction. Int. J. Psychophysiol.; 2015; 95, pp. 101-112. [DOI: https://dx.doi.org/10.1016/j.ijpsycho.2014.04.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24746550]

6. Pires, L.; Leitão, J.; Guerrini, C.; R Simões, M. Event-Related brain potentials in the study of inhibition: Cognitive control, source localization and Age-Related modulations. Neuropsychol. Rev.; 2014; 24, pp. 461-490. [DOI: https://dx.doi.org/10.1007/s11065-014-9275-4]

7. Helmholtz, H. Ueber einige Gesetze der Vertheilung elektrischer Ströme in körperlichen Leitern, mit Anwendung auf die thierisch-elektrischen Versuche (Schluss.). Ann. Phys. Chem.; 1853; 89, pp. 211–233, 353–377. [DOI: https://dx.doi.org/10.1002/andp.18531650603]

8. Papageorgakis, C. Patient Specific Conductivity Models: Characterization of the Skull Bones. Ph.D. Thesis; Université Côte d’Azur: Nice, France, 2017.

9. Vorwerk, J.; Cho, J.; Rampp, S.; Hamer, H.; Knösche, T.; Wolters, C. A guideline for head volume conductor modeling in EEG and MEG. Neuroimage; 2014; 100, pp. 590-607. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2014.06.040]

10. Vorwerk, J.; Aydin, Ü.; Wolters, C.H.; Butson, C.R. Influence of Head Tissue Conductivity Uncertainties on EEG Dipole Reconstruction. Front. Neurosci.; 2019; 13, 531. [DOI: https://dx.doi.org/10.3389/fnins.2019.00531]

11. McCann, H.; Beltrachini, L. Impact of skull sutures, spongiform bone distribution, and aging skull conductivities on the eeg forward and inverse problems. J. Neural Eng.; 2022; 19, 016014. [DOI: https://dx.doi.org/10.1088/1741-2552/ac43f7]

12. Kuratko, D.; Lacik, J.; Koudelka, V.; Vejmola, C.; Wójcik, D.K.; Raida, Z. Forward model of rat electroencephalogram: Comparative study of numerical simulations with measurements on rat head phantoms. IEEE Access; 2022; 10, pp. 92023-92035. [DOI: https://dx.doi.org/10.1109/ACCESS.2022.3202206]

13. Tadel, F.; Baillet, S.; Mosher, J.C.; Pantazis, D.; Leahy, R.M. Brainstorm: A user-friendly application for MEG/EEG analysis. Comput. Intell. Neurosci.; 2011; 2011, 879716. [DOI: https://dx.doi.org/10.1155/2011/879716] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21584256]

14. Oostenveld, R.; Fries, P.; Maris, E.; Schoffelen, J.M. FieldTrip: Open Source Software for Advanced Analysis of MEG, EEG, and Invasive Electrophysiological Data. Comput. Intell. Neurosci.; 2011; 2011, 156869. [DOI: https://dx.doi.org/10.1155/2011/156869] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21253357]

15. Gramfort, A.; Luessi, M.; Larson, E.; Engemann, D.A.; Strohmeier, D.; Brodbeck, C.; Parkkonen, L.; Hämäläinen, M.S. MNE software for processing MEG and EEG data. NeuroImage; 2014; 86, pp. 446-460. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2013.10.027] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24161808]

16. Delorme, A.; Makeig, S. EEGLAB: An open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods; 2004; 134, pp. 9-21. [DOI: https://dx.doi.org/10.1016/j.jneumeth.2003.10.009]

17. Geselowitz, D.B. On Bioelectric Potentials in an Inhomogeneous Volume Conductor. Biophys. J.; 1967; 7, pp. 1-11. [DOI: https://dx.doi.org/10.1016/S0006-3495(67)86571-8]

18. Kybic, J.; Clerc, M.; Abboud, T.; Faugeras, O.; Keriven, R.; Papadopoulo, T. A common formalism for the integral formulations of the forward eeg problem. IEEE Trans. Med. Imaging; 2005; 24, pp. 12-28. [DOI: https://dx.doi.org/10.1109/TMI.2004.837363] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15638183]

19. Ponasso, G.N. A survey on integral equations for bioelectric modeling. Phys. Med. Biol.; 2024; 69, 17TR02. [DOI: https://dx.doi.org/10.1088/1361-6560/ad66a9]

20. Fischl, B.; Dale, A.M. Measuring the thickness of the human cerebral cortex from magnetic resonance images. Proc. Natl. Acad. Sci. USA; 2000; 97, pp. 11050-11055. [DOI: https://dx.doi.org/10.1073/pnas.200033797]

21. Wartman, W.A.; Weise, K.; Rachh, M.; Morales, L.; Deng, Z.; Nummenmaa, A.; Makaroff, S.N. An adaptive h-refinement method for the boundary element fast multipole method for quasi-static electromagnetic modeling. Phys. Med. Biol.; 2024; 69, 055030. [DOI: https://dx.doi.org/10.1088/1361-6560/ad2638]

22. Weise, K.; Wartman, W.A.; Knösche, T.R.; Nummenmaa, A.R.; Makarov, S.N. The effect of meninges on the electric fields in TES and TMS. Numerical modeling with adaptive mesh refinement. Brain Stimul.; 2022; 15, pp. 654-663. [DOI: https://dx.doi.org/10.1016/j.brs.2022.04.009]

23. Partheymüller, P.; Białecki, R.A.; Kuhn, G. Self-adapting algorithm for evaluation of weakly singular integrals arising in the boundary element method. Eng. Anal. Bound. Elem.; 1994; 14, pp. 285-292. [DOI: https://dx.doi.org/10.1016/0955-7997(94)90044-2]

24. Fischl, B. FreeSurfer. NeuroImage; 2012; 62, pp. 774-781. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2012.01.021] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22248573]

25. Nielsen, J.D.; Madsen, K.H.; Puonti, O.; Siebner, H.R.; Bauer, C.; Madsen, C.G.; Saturnino, G.B.; Thielscher, A. Automatic skull segmentation from MR images for realistic volume conductor models of the head: Assessment of the state-of-the-art. NeuroImage; 2018; 174, pp. 587-598. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2018.03.001]

26. Makarov, S.N.; Noetscher, G.M.; Raij, T.; Nummenmaa, A. A Quasi-Static Boundary Element Approach With Fast Multipole Acceleration for High-Resolution Bioelectromagnetic Models. IEEE. Trans. Biomed. Eng.; 2018; 65, pp. 2675-2683. [DOI: https://dx.doi.org/10.1109/TBME.2018.2813261]

27. Miinalainen, T.; Rezaei, A.; Us, D.; Nüßing, A.; Engwer, C.; Wolters, C.H.; Pursiainen, S. A realistic, accurate and fast source modeling approach for the eeg forward problem. NeuroImage; 2019; 184, pp. 56-67. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2018.08.054]

28. Bangera, N.B.; Schomer, D.L.; Dehghani, N.; Ulbert, I.; Cash, S.; Papavasiliou, S.; Eisenberg, S.R.; Dale, A.M.; Halgren, E. Experimental validation of the influence of white matter anisotropy on the intracranial EEG forward solution. J. Comput. Neurosci.; 2010; 29, pp. 371-387. [DOI: https://dx.doi.org/10.1007/s10827-009-0205-z]

29. Güllmar, D.; Haueisen, J.; Reichenbach, J.R. Source analysis of EEG data is an important tool in scientific and clinical applications. This is the first study on EEG source analysis using the BEM. BEM can utilize true point-dipole sources, as opposed to past studies using FEM, which can only use approximate point-dipole sources. NeuroImage; 2010; 51, pp. 145-163.

30. Acar, Z.A.; Makeig, S. Effects of forward model errors on eeg source localization. Brain Topogr.; 2013; 26, pp. 378-396. [DOI: https://dx.doi.org/10.1007/s10548-012-0274-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23355112]

31. Fuchs, M.; Wagner, M.; Kastner, J. Development of volume conductor and source models to localize epileptic foci. J. Clin. Neurophysiol.; 2007; 24, pp. 101-119. [DOI: https://dx.doi.org/10.1097/WNP.0b013e318038fb3e]

32. Whittingstall, K.; Stroink, G.; Gates, L.; Connolly, J.F.; Finley, A. Effects of dipole position, orientation and noise on the accuracy of eeg source localization. Biomed. Eng. Online; 2003; 2, 14. [DOI: https://dx.doi.org/10.1186/1475-925X-2-14]

33. Lew, S.; Wolters, C.H.; Dierkes, T.; Röer, C.; MacLeod, R.S. Accuracy and run-time comparison for different potential approaches and iterative solvers in finite element method based eeg source analysis. Appl. Numer. Math.; 2009; 59, pp. 1970-1988. [DOI: https://dx.doi.org/10.1016/j.apnum.2009.02.006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20161462]

34. Gomez, L.J.; Dannhauer, M.; Koponen, L.M.; Peterchev, A.V. Conditions for numerically accurate tms electric field simulation. Brain Stimul.; 2020; 13, pp. 157-166. [DOI: https://dx.doi.org/10.1016/j.brs.2019.09.015] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31604625]

35. Hämäläinen, M.; Hari, R.; Ilmoniemi, R.J.; Knuutila, J.; Lounasmaa, O.V. Magnetoencephalography—Theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys.; 1993; 65, pp. 413-497. [DOI: https://dx.doi.org/10.1103/RevModPhys.65.413]

36. Schimpf, P.H.; Ramon, C.; Haueisen, J. Dipole models for the eeg and meg. IEEE Trans. Biomed. Eng.; 2002; 49, pp. 409-418. [DOI: https://dx.doi.org/10.1109/10.995679]

37. Buchner, H.; Knoll, G.; Fuchs, M.; Rienäcker, A.; Beckmann, R.; Wagner, M.; Silny, J.; Pesch, J. Inverse localization of electric dipole current sources in finite element models of the human head. Electroencephalogr. Clin. Neurophysiol.; 1997; 102, pp. 267-278. [DOI: https://dx.doi.org/10.1016/S0013-4694(96)95698-9]

38. Olivi, E.; Papadopoulo, T.; Clerc, M. Handling white-matter anisotropy in bem for the eeg forward problem. Proceedings of the 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro; Chicago, IL, USA, 30 March–2 April 2011.

39. Saturnino, G.B.; Puonti, O.; Nielsen, J.D.; Antonenko, D.; Madsen, K.H.; Thielscher, A. SimNIBS 2.1: A Comprehensive Pipeline for Individualized Electric Field Modelling for Transcranial Brain Stimulation. Brain and Human Body Modeling: Computational Human Modeling at EMBC 2018; Springer: Cham, Switzerland, 2019; Chapter 1

40. Van Essen, D.C.; Ugurbil, K.; Auerbach, E.; Barch, D.; Behrens, T.E.J.; Bucholz, R.; Chang, A.; Chen, L.; Corbetta, M.; Curtiss, S.W. et al. The human connectome project: A data acquisition perspective. NeuroImage; 2012; 62, pp. 2222-2231. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2012.02.018]

41. Htet, A.T.; Burnham, E.H.; Noetscher, G.M.; Pham, D.N.; Nummenmaa, A.; Makarov, S.N. Collection of CAD human head models for electromagnetic simulations and their applications. Biomed. Phys. Eng. Express; 2019; 5, 067005. [DOI: https://dx.doi.org/10.1088/2057-1976/ab4c76]

42. Gabriel, C. Compilation of the Dielectric Properties of Body Tissues at Rf and Microwave Frequencies. 1996; N.AL/OE-TR-1996-0037, Occupational and Environmental Health Directorate, Radiofrequency Radiation Division, Brooks Air Force Base, Texas (USA) Available online: http://niremf.ifac.cnr.it/docs/DIELECTRIC/home.html (accessed on 1 September 2024).

43. Hasgall, P.A.; Gennaro, F.D.; Baumgartner, C.; Neufeld, E.; Lloyd, B.; Gosselin, M.C.; Payne, D.; Klingenböck, A.; Kuster, N. IT’IS Database for Thermal and Electromagnetic Parameters of Biological Tissues, 2 2022. Version 4.1. Available online: https://itis.swiss/virtual-population/tissue-properties/database/low-frequency-conductivity/ (accessed on 1 September 2024).

44. Vorwerk, J.; Wolters, C.H.; Baumgarten, D. Global sensitivity of EEG source analysis to tissue conductivity uncertainties. Front. Hum. Neurosci.; 2024; 18, 1335212. [DOI: https://dx.doi.org/10.3389/fnhum.2024.1335212]

45. Akhtari, M.; Bryant, H.C.; Mamelak, A.N.; Flynn, E.R.; Heller, L.; Shih, J.J.; Mandelkem, M.; Matlachov, A.; Ranken, D.M.; Best, E.D. et al. Conductivities of three-layer live human skull. Brain Topogr.; 2002; 14, pp. 151-167. [DOI: https://dx.doi.org/10.1023/A:1014590923185]

46. Baumann, S.B.; Wozny, D.R.; Kelly, S.K.; Meno, F.M. The electrical conductivity of human cerebrospinal fluid at body temperature. IEEE Trans. Biomed. Eng.; 1997; 44, pp. 220-223. [DOI: https://dx.doi.org/10.1109/10.554770]

47. Dannhauer, M.; Lanfer, B.; Wolters, C.H.; Knösche, T.R. Modeling of the human skull in EEG source analysis. Hum. Brain Mapp.; 2011; 32, pp. 1383-1399. [DOI: https://dx.doi.org/10.1002/hbm.21114] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20690140]

48. Haueisen, J.; Ramon, C.; Eiselt, M.; Brauer, H.; Nowak, H. Influence of tissue resistivities on neuromagnetic fields and electric potentials studied with a finite element model of the head. IEEE Trans. Biomed. Eng.; 1997; 44, pp. 727-735. [DOI: https://dx.doi.org/10.1109/10.605429]

49. Ramon, C.; Schimpf, P.; Haueisen, J.; Holmes, M.; Ishimaru, A. Role of soft bone, CSF and gray matter in EEG simulations. Brain Topogr.; 2004; 16, pp. 245-248. [DOI: https://dx.doi.org/10.1023/B:BRAT.0000032859.68959.76]

50. Hirata, A.; Niitsu, M.; Phang, C.R.; Kodera, S.; Kida, T.; Rashed, E.A.; Fukunaga, M.; Sadato, N.; Wasaka, T. High-resolution eeg source localization in personalized segmentation-free head model with multi-dipole fitting. Phys. Med. Biol.; 2024; 69, 055013. [DOI: https://dx.doi.org/10.1088/1361-6560/ad25c3]

51. Opitz, A.; Paulus, W.; Will, S.; Antunes, A.; Thielscher, A. Determinants of the electric field during transcranial direct current stimulation. NeuroImage; 2015; 109, pp. 140-150. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2015.01.033]

52. Wagner, T.A.; Zahn, M.; Grodzinsky, A.J.; Pascual-Leone, A. Three-dimensional head model simulation of transcranial magnetic stimulation. IEEE Transact. Biomed. Eng.; 2004; 51, pp. 1586-1598. [DOI: https://dx.doi.org/10.1109/TBME.2004.827925]

53. Cignoni, P.; Callieri, M.; Corsini, M.; Dellepiane, M.; Ganovelli, F.; Ranzuglia, G. MeshLab: An Open-Source Mesh Processing Tool. Eurographics Italian Chapter Conference, Salerno, Italy; Scarano, V.; Chiara, R.D.; Erra, U. The Eurographics Association: Crete, Greece, 2008.

54. Kazhdan, M.; Hoppe, H. Screened Poisson surface reconstruction. ACM Trans. Graph. (TOG); 2013; 32, 29. [DOI: https://dx.doi.org/10.1145/2487228.2487237]

55. Taubin, G. Curve and surface smoothing without shrinkage. Proceedings of the IEEE International Conference on Computer Vision; Cambridge, MA, USA, 20–23 June 1995; pp. 852-857.

56. Pittau, F.; Grouiller, F.; Spinelli, L.; Seeck, M.; Michel, C.M.; Vulliemoz, S. The Role of Functional Neuroimaging in Pre-Surgical Epilepsy Evaluation. Front. Neurol.; 2014; 5, 31. [DOI: https://dx.doi.org/10.3389/fneur.2014.00031] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24715886]

57. Ebersole, J.S.; Hawes-Ebersole, S. Clinical Application of Dipole Models in the Localization of Epileptiform Activity. J. Clin. Neurophysiol.; 2007; 24, pp. 120-129. [DOI: https://dx.doi.org/10.1097/WNP.0b013e31803ece13]

58. Fiedler, P.; Fonseca, C.; Supriyanto, E.; Zanow, F.; Haueisen, J. A high-density 256-channel cap for dry electroencephalography. Human Brain Mapp.; 2022; 43, pp. 1295-1308. [DOI: https://dx.doi.org/10.1002/hbm.25721]

59. Graichen, U.; Eichardt, R.; Fiedler, P.; Strohmeier, D.; Zanow, F.; Haueisen, J. SPHARA—A Generalized Spatial Fourier Analysis for Multi-Sensor Systems with Non-Uniformly Arranged Sensors: Application to EEG. PLoS ONE; 2015; 10, e0121741. [DOI: https://dx.doi.org/10.1371/journal.pone.0121741] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25885290]

60. Barnard, A.C.; Duck, I.M.; Lynn, M.S. The application of electromagnetic theory to electrocardiology. i. derivation of the integral equations. Biophys. J.; 1967; 7, pp. 443-462. [DOI: https://dx.doi.org/10.1016/S0006-3495(67)86598-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/6048873]

61. Gelernter, H.L.; Swihart, J.C. A Mathematical-Physical Model of the Genesis of the Electrocardiogram. Biophys. J.; 1964; 4, pp. 285-301. [DOI: https://dx.doi.org/10.1016/S0006-3495(64)86783-7]

62. Greengard, L.; Rokhlin, J. A fast algorithm for particle simulations. J. Comp. Phys.; 1987; 73, pp. 325-348. [DOI: https://dx.doi.org/10.1016/0021-9991(87)90140-9]

63. Nunez, P.L.; Srinivasan, R. Electric Fields of the Brain: The Neurophysics of EEG; 2nd ed. Oxford University Press: Oxford, UK, 2006.

64. Makarov, S.N.; Wartman, W.A.; Noetscher, G.M.; Fujimoto, K.; Zaidi, T.; Burnham, E.T.; Daneshzand, M.; Nummenmaa, A. Degree of improving TMS focality through a geometrically stable solution of an inverse TMS problem. NeuroImage; 2021; 241, 118437. [DOI: https://dx.doi.org/10.1016/j.neuroimage.2021.118437] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34332043]

65. Makarov, S.N.; Hämäläinen, M.; Yoshio, O.; Noetscher, G.M.; Jyrki, A.; Nummenmaa, A. Boundary element fast multipole method for enhanced modeling of neurophysiological recordings. IEEE Trans. Biomed. Eng.; 2021; 68, pp. 308-318. [DOI: https://dx.doi.org/10.1109/TBME.2020.2999271]

66. Makarov, S.N.; Noetscher, G.M.; Nazarian, A. Low-Frequency Electromagnetic Modeling for Electrical and Biological Systems Using MATLAB; Wiley: New York, NY, USA, 2015.

67. Rick, B.; Leslie, G. A Short Course on Fast Multipole Methods. 1997; Available online: https://math.nyu.edu/~greengar/shortcourse_fmm.pdf (accessed on 1 September 2024).

68. Rokhlin, V. Rapid solution of integral equations of classical potential theory. J. Comput. Phys.; 1985; 60, pp. 187-207. [DOI: https://dx.doi.org/10.1016/0021-9991(85)90002-6]

69. Askham, T.; Gimbutas, Z.; Greengard, L.; Lu, L.; Magland, J.; Malhotra, D.; O’Neil, M.; Rachh, M.; Rokhlin, V.; Vico, F. FMM3D: A Fast Multipole Method Library for Three-Dimensional Problems. Available online: https://github.com/flatironinstitute/FMM3D (accessed on 1 September 2024).

70. Feischl, M.; Führer, T.; Heuer, N.; Karkulik, M.; Praetorius, D. Adaptive Boundary Element Methods. Arch. Computat. Methods Eng.; 2015; 22, pp. 309-389. [DOI: https://dx.doi.org/10.1007/s11831-014-9114-z]

71. Wartman, W.A.; Ponasso, G.N.N.; Qi, Z.; Haueisen, J.; Maess, B.; Knösche, T.R.; Weise, K.; Noetscher, G.M.; Raij, T.; Makaroff, S.N. Fast and accurate eeg/meg bem-based forward problem solution for high-resolution head models. bioRxiv; 2024; [DOI: https://dx.doi.org/10.1101/2024.06.07.598024]

72. Scherg, M. Fundamentals of dipole source potential analysis. Adv. Audiol.; 1990; 6, 25.

73. Phillips, C. Source Localisation in EEG: Combining Anatomical and Functional Constraints. Ph.D. Thesis; Université de Liège: Liège, Belgium, 2001.

74. Sarvas, J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys. Med. Biol.; 1987; 32, pp. 11-22. [DOI: https://dx.doi.org/10.1088/0031-9155/32/1/004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3823129]

75. Piastra, M.C.; Oostenveld, R.; Homölle, S.; Han, B.; Chen, Q.; Oostendorp, T. How to assess the accuracy of volume conduction models? A validation study with stereotactic eeg data. Front. Hum. Neurosci.; 2024; 18, 1279183. [DOI: https://dx.doi.org/10.3389/fnhum.2024.1279183] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38410258]