Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2791331

Formats

Article sections

Authors

Related links

Appl Numer Math. Author manuscript; available in PMC 2010 August 1.

Published in final edited form as:

Appl Numer Math. 2009 August; 59(8): 1970–1988.

doi: 10.1016/j.apnum.2009.02.006PMCID: PMC2791331

NIHMSID: NIHMS120338

See other articles in PMC that cite the published article.

Accuracy and run-time play an important role in medical diagnostics and research as well as in the field of neuroscience. In Electroencephalography (EEG) source reconstruction, a current distribution in the human brain is reconstructed noninvasively from measured potentials at the head surface (the EEG inverse problem). Numerical modeling techniques are used to simulate head surface potentials for dipolar current sources in the human cortex, the so-called EEG forward problem.

In this paper, the efficiency of algebraic multigrid (AMG), incomplete Cholesky (IC) and Jacobi preconditioners for the conjugate gradient (CG) method are compared for iteratively solving the finite element (FE) method based EEG forward problem. The interplay of the three solvers with a full subtraction approach and two direct potential approaches, the Venant and the partial integration method for the treatment of the dipole singularity is examined. The examination is performed in a four-compartment sphere model with anisotropic skull layer, where quasi-analytical solutions allow for an exact quantification of computational speed versus numerical error. Specifically-tuned constrained Delaunay tetrahedralization (CDT) FE meshes lead to high accuracies for both the full subtraction and the direct potential approaches. Best accuracies are achieved by the full subtraction approach if the homogeneity condition is fulfilled. It is shown that the AMG-CG achieves an order of magnitude higher computational speed than the CG with the standard preconditioners with an increasing gain factor when decreasing mesh size. Our results should broaden the application of accurate and fast high-resolution FE volume conductor modeling in source analysis routine.

Electroencephalography (EEG) based source reconstruction of cerebral activity (the EEG inverse problem) is an important tool both in clinical practice and research [35,23], and in cognitive neuroscience [2]. Methods for solving the inverse problem are based on solutions to the corresponding forward problem, i.e., the simulation of EEG potentials for a given primary source in the brain using a volume-conduction model of the human head. While the theory of this forward problem is well established and many numerical implementations exist, there remain unresolved questions regarding the accuracy and efficiency of contemporary approaches. In this study, we compared a range of numerical techniques and source representation approaches and have shown that careful choice of both are critical in order to solve realistic electroencephalographic forward (and inverse) problems.

The general approach for solving bioelectric field problems under realistic conditions is well established. All quantitative solutions for the EEG forward problem are based on the quasi-static Maxwell equations [25]. The primary sources are electrolytic currents within the dendrites of the large pyramidal cells of activated neurons in the human cortex. Even if there are also smoother models [33], most often the primary sources are formulated as a mathematical point current dipole [25,6,18]. The finite element (FE) method is often used for the solution of the forward problem, because it allows for a realistic representation of the complicated head volume conductor with its tissue conductivity inhomogeneities and anisotropies [45,3,1,34,4,15,19,37,26,39,7].

To implement the point current dipole as a current source in the brain, the FE method requires careful consideration of the singularity of the potential at the source position. One way to address the singularity is to use a *subtraction approach*, which divides the total potential into an analytically known singularity potential and a singularity-free correction potential, which can then be approximated numerically using an FE approach [3,1,34,15,26,41,39]. For the correction potential, the existence and uniqueness for a weak solution in a zero-mean function space have been proven and FE convergence properties are known [39]. The subtraction FE approach has thus a sound mathematical basis for point current dipole models. Another family of source representation methods, known as *direct FE approaches* to the total potential [45,1,4,37,26], are computationally less expensive, but also mathematically less sound under the assumption that a point dipole is the more realistic source model. In our previous work [41,44], we compared the *projected subtraction approach* [39] with the two direct approaches using *partial integration* [45,1,37] and *Venant* [4] in regular and geometry-adapted 2mm hexahedral FE meshes of a multi-layer sphere model for which quasi-analytical solutions exist [5]. The Venant approach was found to be the overall best choice in FE source analysis practice [41,44]. It was however speculated that an improved numerical quadrature might solve the accuracy problems of the projected subtraction approach for eccentric sources, so that the highest future potential was given to an improved implementation of the FE subtraction approach [44]. It has recently been shown that a *full subtraction approach* [7] using an improved numerical quadrature leads to an order of magnitude more accurate solution than the projected subtraction approach [39], especially when considering sources that are close to a conductivity inhomogeneity.

Another general prerequisite for FE modeling of bioelectric fields is the generation of a mesh that represents the geometry and electric properties of the volume conductor. An effective meshing strategy will balance acceptable forward problem accuracy against reasonable computation times and memory usage. Very high accuracies can be achieved by making use of a Constrained Delaunay Tetrahedralization (CDT) [30,29] in combination with a full subtraction approach [7]. Adaptive methods, using local refinement around the source singularity [3,34], are another potential utility but they preclude the use of fast transfer matrices [37,8,43,12] and lose efficiency in solving the inverse problem (see discussion section).

Solving the forward problem is rarely the ultimate goal in calculating bioelectric fields but rather a step towards solving the associated inverse problem. Thus the quest for numerical accuracy and efficiency of the forward solution requires some anticipation of the ultimate use in inverse solutions. The longtime state-of-the-art approach has been to solve an FE equation system for each anatomically and physiologically meaningful dipolar source (each source results in one FE right-hand side (RHS) vector) [3,1,34,4,15]. The use of standard direct (banded LU factorization for a 2D source analysis scenario [1]) or iterative (Conjugate Gradient (CG) without preconditioning [3] or Successive OverRelaxation (SOR) [26]) FE solver techniques limit the overall resolution of the geometric model because of their computational cost. The preconditioned CG method was used with standard preconditioners like Jacobi (Jacobi-CG) [37] or incomplete Cholesky without fill-in, IC(0)-CG [4].

One recent approach to achieve efficient computation of the FE-based forward problem is to pre-compute transfer matrices that encapsulate the relationship between source locations and sensor sites based only on the geometric and conductivity characteristics of the volume conductor, i.e., they are independent of the source. Techniques exist to construct transfer matrices for problem formulations based on EEG [37,12] or combined EEG and MEG [8,43]. Using this principle, for each head model, one only has to solve one large sparse FE system of equations for each of the possible sensor locations in order to compute the full transfer matrix. Each forward solution is then reduced to multiplication of the transfer matrix by an FE RHS vector containing the source load. Exploiting the fact that the number of sensors (currently up to about 600) is much smaller than the number of reasonable dipolar sources (tens of thousands), the transfer matrix approach is substantially faster than the state-of-the-art forward approach (i.e., solving an FE equation system for each source) and can be applied to inverse reconstruction algorithms in both continuous and discrete source parameter space for EEG and MEG. Still, the solution of hundreds of large linear FE equation systems for the construction of the transfer matrices is a major time consuming part within FE-based source analysis.

The first goal of this study was therefore to compare the numerical accuracy of the full subtraction approach [7] with the two direct approaches using partial integration [45,1,37] and Venant [4] in specifically-tuned CDT meshes of an anisotropic four-compartment sphere model. We then examine the interplay of the source model approaches with three FE solver methods: a Jacobi-CG, an incomplete Cholesky CG (e.g., [27]), and an algebraic multigrid preconditioned CG (AMG-CG), which has already shown to be especially suited for problems with discontinuous and anisotropic coefficients [22,32,21,9,38].

In the quasi-static approximation of the Maxwell equations, the distribution of electric potentials Φ in the head domain Ω of conductivity *σ*, resulting from a primary current **j**^{p} is governed by the Poisson equation with homogeneous Neumann boundary conditions on the head surface Γ = Ω [20,25], which we can express as

$$\nabla \cdot (\sigma \nabla \Phi )=\nabla \cdot {\mathbf{j}}^{p}={\mathrm{J}}^{p}\phantom{\rule{thickmathspace}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\Omega ,\phantom{\rule{1em}{0ex}}\langle \sigma \nabla \Phi ,\mathbf{n}\rangle =0\phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma ,$$

(1)

with **n** the unit surface normal, and assuming a reference electrode with given potential, i.e., Φ(**x**_{ref}) = 0. The primary currents are modeled by a mathematical dipole at position ${\mathbf{x}}_{0}\in {\mathbb{R}}^{3}$ with moment ${\mathbf{M}}_{0}\in {\mathbb{R}}^{3}$ [25,6,18],

$${\mathrm{J}}^{p}=\nabla \cdot {\mathbf{j}}^{p}\left(\mathbf{x}\right)=\nabla \cdot \left({\mathbf{M}}_{0}\delta (\mathbf{x}-{\mathbf{x}}_{0})\right),$$

(2)

where *δ* is the Dirac delta distribution.

One of the key questions for all three-dimensional EEG forward modeling techniques is the appropriate treatment of the potential singularity introduced into the differential equation by the formulation of the mathematical dipole (2). This study examined the interplay of FE solver methods (see Section 2.2) with the solution accuracy in four-layer sphere models applying three singularity treatment techniques: a full subtraction approach, a partial integration direct method and a Venant direct method.

The subtraction approach [3,1,39,7] splits the total potential Φ into two parts,

$$\Phi ={\Phi}_{0}+{\Phi}_{\mathrm{corr}},$$

(3)

where the singularity potential, Φ_{0}, is defined as the solution for a dipole in an unbounded homogeneous conductor with constant conductivity σ_{0}. ${\sigma}_{0}\in {\mathbb{R}}^{3\times 3}$ is the conductivity at the source position, which is assumed to be constant in a non-empty subdomain Ω_{0} around **x**_{0}, in the following called the *homogeneity condition*. The solution of Poisson’s equation under these conditions for the singularity potential

$$\nabla \cdot ({\sigma}_{0}\nabla {\Phi}_{0})=\nabla \cdot {\mathbf{j}}^{p}$$

(4)

can be formed analytically for the mathematical dipole (2) [7] as

$${\Phi}_{0}\left(\mathbf{x}\right)=\frac{1}{4\pi \sqrt{\mathrm{det}\phantom{\rule{thinmathspace}{0ex}}{\sigma}_{0}}}\frac{\langle {\mathbf{M}}_{0},{\left({\sigma}_{0}\right)}^{-1}(\mathbf{x}-{\mathbf{x}}_{0})\rangle}{{\langle {\left({\sigma}_{0}\right)}^{-1}(\mathbf{x}-{\mathbf{x}}_{0}),(\mathbf{x}-{\mathbf{x}}_{0})\rangle}^{3\u22152}}.$$

(5)

Subtracting (4) from (1) yields a Poisson equation for the correction potential

$$-\nabla \cdot (\sigma \nabla {\Phi}_{\mathrm{corr}})=-\nabla \cdot (({\sigma}_{0}-\sigma )\nabla {\Phi}_{0})\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{thickmathspace}{0ex}}\Omega ,$$

(6)

with inhomogeneous Neumann boundary conditions at the surface:

$$\langle \sigma \nabla {\Phi}_{\mathrm{corr}},\mathbf{n}\rangle =-\langle \sigma \nabla {\Phi}_{0},\mathbf{n}\rangle \phantom{\rule{1em}{0ex}}\text{on}\phantom{\rule{thickmathspace}{0ex}}\Gamma .$$

(7)

The advantage of (6) is that the right-hand side is free of any source singularity, because of the homogeneity condition — the conductivity *σ*_{0} – *σ* is zero in Ω_{0}. Existence and uniqueness of the solution and FE convergence properties are shown for the correction potential in [39]. For the numerical approximation of the correction potential, we use the FE method with piecewise linear basis functions * _{i}*. When projecting the correction potential into the FE space, i.e., ${\Phi}_{\mathrm{corr}}\left(\mathbf{x}\right)\approx {\Phi}_{\mathrm{corr},h}\left(\mathbf{x}\right)={\sum}_{j=1}^{{N}_{h}}{\phi}_{j}\left(\mathbf{x}\right){\underset{\u2012}{u}}_{\mathrm{corr},h}^{\left[j\right]}$ (

$${K}_{h}{\underset{\u2012}{u}}_{\mathrm{corr},h}={\underset{\u2012}{j}}_{\mathrm{corr},h},$$

(8)

with the stiffness matrix

$${K}_{h}^{[i,j]}={\int}_{\Omega}\langle \sigma \nabla {\phi}_{j},\nabla {\phi}_{i}\rangle dx,$$

(9)

for ${K}_{h}\in {\mathbb{R}}^{{N}_{h}\times {N}_{h}}$, and the right-hand side vector ${\underset{\u2012}{j}}_{\mathrm{corr},h}\in {\mathbb{R}}^{{N}_{h}}$ with entries

$${\underset{\u2012}{j}}_{\mathrm{corr},h}^{\left[i\right]}={\int}_{\Omega}\langle ({\sigma}_{0}-\sigma )\nabla {\Phi}_{0},\nabla {\phi}_{i}\left(x\right)\rangle dx-{\int}_{\partial \Omega}{\phi}_{i}\left(x\right)\langle n\left(x\right),{\sigma}_{0}\nabla {\Phi}_{0}\left(x\right)\rangle dx.$$

(10)

We then seek for the coefficient vector ${\underset{\u2012}{u}}_{\mathrm{corr},h}=({\underset{\u2012}{u}}_{\mathrm{corr},h}^{\left[1\right]},\dots ,{\underset{\u2012}{u}}_{\mathrm{corr},h}^{\left[{N}_{h}\right]})\in {\mathbb{R}}^{{N}_{h}}$ and, using (3), compute the total potential. In [7], the theoretical reasoning and a validation in a four-compartment sphere model with anisotropic skull is given for the fact that second order integration is necessary and sufficient for the right-hand side integration in Equation (10). Direct comparisons with the projected subtraction approach from [39] have shown that the full subtraction approach is an order of magnitude more accurate for dipole sources close to a conductivity discontinuity [7].

Multiplying both sides of Equation (1) by a linear FE basis function * _{i}* and integrating over the head domain leads to a partial integration direct approach for the total potential [1,36,17] expressed as

$${\int}_{\Omega}\nabla \cdot (\sigma \nabla \Phi ){\phi}_{i}dx={\int}_{\Omega}\nabla \cdot {\mathbf{j}}^{p}{\phi}_{i}dx.$$

Integration by parts, applied to both sides of the above equation, yields

$$-{\int}_{\Omega}\langle \sigma \nabla \Phi ,\nabla {\phi}_{i}\rangle dx+{\int}_{\Gamma}\langle \sigma \nabla \Phi ,\mathbf{n}\rangle {\phi}_{i}d\Gamma =-{\int}_{\Omega}\langle {\mathbf{j}}^{p},\nabla {\phi}_{i}\rangle dx+{\int}_{\Gamma}\langle {\mathbf{j}}^{p},\mathbf{n}\rangle {\phi}_{i}d\Gamma .$$

Using the homogeneous Neumann boundary condition from Equation (1) and the fact that the current density vanishes on the head surface, we arrive at

$${\int}_{\Omega}\langle \sigma \nabla \Phi ,\nabla {\phi}_{i}\rangle dx={\int}_{\Omega}\langle {\mathbf{j}}^{p},\nabla {\phi}_{i}\rangle dx\stackrel{\left(2\right)}{=}\langle {\mathbf{M}}_{0},\nabla {\phi}_{i}\left({\mathbf{x}}_{0}\right)\rangle .$$

Setting $\Phi \left(\mathbf{x}\right)\approx {\Phi}_{h}\left(\mathbf{x}\right)={\sum}_{j=1}^{{N}_{h}}{\phi}_{j}\left(\mathbf{x}\right){\underset{\u2012}{u}}_{h}^{\left[j\right]}$, leads to the linear system

$${K}_{h}{\underset{\u2012}{u}}_{h}={\underset{\u2012}{j}}_{\mathrm{PI},h},$$

(11)

with the same stiffness matrix as in (9) and the right-hand side vector ${\underset{\u2012}{j}}_{\mathrm{PI},h}\in {\mathbb{R}}^{{N}_{h}}$ with entries

$${\underset{\u2012}{j}}_{\mathrm{PI},h}^{\left[i\right]}=\{\begin{array}{cc}\hfill \langle {\mathbf{M}}_{0},\nabla {\phi}_{i}\left({\mathbf{x}}_{0}\right)\rangle \hfill & \hfill \text{if}\phantom{\rule{thickmathspace}{0ex}}i\in \text{NODESOFELE}\left({\mathbf{x}}_{0}\right),\hfill \\ \hfill 0\hfill & \hfill \text{otherwise}.\hfill \end{array}\phantom{\}}$$

(12)

The function NODESOFELE(**x**_{0}) determines the set of nodes of the element which contains the dipole at position **x**_{0}. Note that while the right-hand side vector (10) is fully populated, *j*_{PI,h} has only |NODESOFELE| non-zero entries. Here, |·| denotes the number of elements in the set NODESOFELE. For the linear basis functions * _{i}* considered here, the right-hand side (12) and thus the computed solution for the total potential in (11) will be constant for all

The Venant potential approach [4] follows the principle of Saint Venant and is made up from monopolar loads on all neighboring FE nodes so that the dipolar moment is fulfilled and the source load is as regular as possible. In this case, variational and FE techniques yield the linear system

$${K}_{h}{\underset{\u2012}{u}}_{h}={\underset{\u2012}{j}}_{\text{Venant},h}$$

(13)

with the same stiffness matrix as in (9). The right-hand side vector ${\underset{\u2012}{j}}_{\text{Venant},h}\in {\mathbb{R}}_{h}^{N}$ has only *C* non-zero entries, if *C* is the number of neighboring FE nodes to that FE node which is closest to the location of the dipole [4].

The solution of hundreds of large scale systems of equations (8), (11) or (13) with the same symmetric positive definite (SPD) stiffness matrix (9) is the major time consuming task of the inverse source localization process. The spectral condition of the SPD matrix *K _{h}* is equal to

$${\kappa}_{2}\left({K}_{h}\right)=\frac{{\lambda}_{\mathrm{max}}}{{\lambda}_{\mathrm{min}}}$$

with λ_{max} the largest and λ_{min} the smallest eigenvalues, respectively, of *K _{h}* [11, §2.10]. The condition number behaves asymptotically as $\mathcal{O}\left({h}^{-2}\right)$ and condition numbers of more than 10

The Preconditioned Conjugate Gradient (PCG) iterative solver shown in Algorithm 1 (see, e.g.,[27,11,24]) can provide efficient procedures for such problems. Note that, in theory, the convergence speed of the PCG is independent of the right-hand side *j _{h}* of the linear equation system [11, §3.4]. The goal of a preconditioner, ${C}_{h}\in {\mathbb{R}}^{{N}_{h}\times {N}_{h}}$, is the reduction of ${\kappa}_{2}\left({C}_{h}^{-1}{K}_{h}\right)$ for the preconditioned equation system ${C}_{h}^{-1}{K}_{h}{\underset{\u2012}{u}}_{h}={C}_{h}^{-1}{\underset{\u2012}{j}}_{h}$. Further requirements are that it is cheap with regard to arithmetic and memory costs to solve linear systems

**Theorem 2.1 (Error estimate for PCG method)** Let K_{h} and C_{h} be positive definite. If ${\underset{\u2012}{u}}_{h}^{\ast}$ denotes the exact solution of the equation system, then the k^{th} iterate of the PCG method ${\underset{\u2012}{u}}_{h}^{k}$ fulfills the following energy norm estimate

$${\Vert {\underset{\u2012}{u}}_{h}^{k}-{\underset{\u2012}{u}}_{h}^{\ast}\Vert}_{{K}_{h}}\le {c}^{k}=\frac{2}{1+{c}^{2k}}{\Vert {\underset{\u2012}{u}}_{h}^{0}-{\underset{\u2012}{u}}_{h}^{\ast}\Vert}_{{K}_{h}},\phantom{\rule{1em}{0ex}}c\u2254\frac{\sqrt{{\kappa}_{2}\left({C}_{h}^{-1}{K}_{h}\right)}-1}{\sqrt{{\kappa}_{2}\left({C}_{h}^{-1}{K}_{h}\right)+1}}.$$

**Proof:** Hackbusch [11, Theorem 9.4.14].

As indicated in Algorithm 1, the PCG method is stopped after the *k ^{th}* iteration if the relative error, i.e., ${\underset{\u2012}{e}}_{h}^{k}={\underset{\u2012}{u}}_{h}^{k}-{\underset{\u2012}{u}}_{h}^{\ast}$ in the controllable ${K}_{h}{C}_{h}^{-1}{K}_{h}\text{-energy}$ norm is below a given ACCURACY. Besides the scalar products, the sum of vectors and the multiplication of a vector with a scalar value, the most important steps in Algorithm 1 are the multiplication of the sparse stiffness matrix

In the following, the three different preconditioners for the CG method are presented and their relative performances evaluated in Section 4.

The simplest preconditioner is the scaling or Jacobi-preconditioning ([24, pp.265f], [27, pp.257f]), where

$${C}_{h}\u2254{D}_{h}^{2},\phantom{\rule{1em}{0ex}}{D}_{h}\u2254\mathrm{DIAG}(\sqrt{{K}_{h}^{\left[11\right]}},\dots ,\sqrt{{K}_{h}^{\left[{N}_{h}{N}_{h}\right]}}).$$

When splitting the Jacobi-preconditioner between left and right (row and column scaling), one has to solve ${\stackrel{~}{K}}_{h}{\underset{\u2012}{v}}_{h}={D}_{h}^{-1}{\underset{\u2012}{j}}_{h}$ with ${\stackrel{~}{K}}_{h}={D}_{h}^{-1}{K}_{h}{D}_{h}^{-\mathit{tr}}$ and ${\underset{\u2012}{u}}_{h}={D}_{h}^{-\mathit{tr}}{\underset{\u2012}{v}}_{h}$. Row and column scaling preserves symmetry, so that the scaled matrix ${\stackrel{~}{K}}_{h}$ is again SPD with unit diagonal entries. The scaling may lead to a first substantial condition improvement, because diagonal entries in *K _{h}* of FE nodes from inside the skull are much smaller than from outside (because of a jump in conductivity at each internal and external boundary) and because it can be shown, that the smallest (largest) eigenvalue of a symmetric matrix is at most (at least) as large as the smallest (largest) diagonal element, so that the condition number is at least as large as the quotient of maximal and minimal diagonal element [27, p.258].

**Theorem 2.2** Let K_{h} be SPD and ${C}_{h}\u2254{D}_{h}^{2}$ the Jacobi-preconditioner. Assume that each row of K_{h} does not contain more than d nonzero entries. Then, for all diagonal matrices ${\stackrel{~}{D}}_{h}^{-1}$, it is

$${\kappa}_{2}\left({C}_{h}^{-1}{K}_{h}\right)\le \mathrm{d}\phantom{\rule{thickmathspace}{0ex}}{\kappa}_{2}\left({\stackrel{~}{D}}_{h}^{-1}{K}_{h}\right),$$

i.e., the chosen diagonal preconditioner is close to the optimal one.

**Proof:** Hackbusch [11, Theorem 8.3.3].

The SPD stiffness matrix *K _{h}* can be decomposed into a left triangular matrix

$${\stackrel{\u02d8}{K}}_{h}={\mathrm{Id}}_{h}+\frac{1}{1+\varsigma}({E}_{h}+{E}_{h}^{\mathit{tr}}).$$

(14)

For sufficiently large $\varsigma \in {\mathbb{R}}_{0}^{+}$, the existence of IC(0) is guaranteed, but with increasing ς, the preconditioning effect decreases. Note that for certain special cases, a condition improvement to $\mathcal{O}\left({h}^{-1}\right)$ can be proven as, e.g., when using a modified ILU_{ω}-preconditioning with *ω* = −1 (in the symmetric case, the ILU_{0} is equal to the IC(0)) for diagonally dominant symmetric matrices arising from a 5-point discretization of a two-dimensional Poisson equation (Hackbusch [11, Theorem 8.5.15 and Remarks 8.5.16,17]).

The above preconditioning methods have the disadvantage that the convergence rate, i.e., the factor by which the error is reduced in each iteration, is still dependent on the mesh size *h*. With decreasing mesh size and thus increasing order of the equation system, the convergence rate tends to 1 from below, so that the number of iterations needed to achieve a given accuracy increases. For the Geometric Multi-Grid (GMG), an *h*-independent convergence rate *ρ* < 1 and an *h*-independent condition number has been proven in [11, Lemma 10.7.1,Theorem 10.7.15] as

$${\kappa}_{2}\left({C}_{h}^{-1}{K}_{h}\right)\le \frac{1}{1-{\rho}^{m}},$$

(15)

with *C _{h}* the preconditioner resulting from

In contrast to GMG, in which a grid hierarchy is required explicitly, Algebraic MG (AMG) is able to construct a matrix hierarchy and corresponding transfer operators based only on the entries in *K _{h}* (see, e.g., [22,32,21,9]). It is well known that the classical AMG method is robust for M-matrices and, with regard to our application, that small positive off-diagonal entries are admissible [22,32,21]. The method is especially well suited for our problem with discontinuous and anisotropic coefficients, in which an optimal tuning of the GMG is difficult ([22, §4.1,4.6.4],[32, §4.1]). Stand-alone AMG is hardly ever optimal as there may be some very specific error components which are reduced with significantly less efficiency, causing a few eigenvalues of the AMG iteration matrix to be much closer to 1 than the remaining ones [32, §3.3]. In such a case, acceleration by means of using AMG as a basis for the CG method eliminates these particular frequencies very efficiently.

As in GMG, the basic idea in AMG is to reduce high and low frequency components of the error by the efficient interplay of smoothing and coarse grid correction, respectively. In analogy to GMG, the denotation *coarse grids* will be used, although these are purely virtual and do not have to be constructed explicitly as FE meshes. The diagonal entry of the *i ^{th}* row of

- Coarsening: define the splitting
*ω*=_{h}*ω*_{C}*ω*of_{F}*ω*into sets of coarse and fine grid nodes_{h}*ω*and_{C}*ω*, respectively._{F} - Transfer operators: prolongation ${P}_{h,H}:{\mathbb{R}}^{{N}_{H}}\mapsto {\mathbb{R}}^{{N}_{h}}$ and its adjoint as the restriction$${R}_{H,h}\u2254{P}_{h,H}^{\mathit{tr}}.$$(16)
- Definition of the coarse matrix by Galerkin’s method, i.e.,Because of (b), ${K}_{H}\in {\mathbb{R}}^{{N}_{H}\times {N}_{H}}$ is again SPD.$${K}_{H}\u2254{R}_{H,h}{K}_{h}{P}_{h,H}.$$(17)

The coarsening process has the task of reducing the number of nodes such that *N _{H}* = |

$${\mathit{Ne}}_{h}^{i}=\left\{j\mid \phantom{\rule{thinmathspace}{0ex}}\mid {K}_{h}^{\left[ij\right]}\mid \ge \zeta \mid {K}_{h}^{[i,i]}\mid ,i\ne j\right\},$$

(18)

$$\begin{array}{cc}\hfill {S}_{h}^{i}& =\left\{j\in {\mathit{Ne}}_{h}^{i}\mid \phantom{\rule{thinmathspace}{0ex}}\mid {K}_{h}^{\left[ij\right]}\mid >\text{coarse}(i,j,{K}_{h})\right\},\hfill \\ \hfill {S}_{h}^{i,T}& =\left\{j\in {\mathit{Ne}}_{h}^{i}\mid i\in {S}_{h}^{j})\right\},\hfill \end{array}$$

(19)

where $N{e}_{h}^{i}$ is the index set of neighbors (a pre-selection is carried out by the threshold-parameter $\varsigma \in {\mathbb{R}}_{0}^{+}$), ${S}_{h}^{i}$ denotes the index set of nodes with a *strong connection* from node *i* and ${S}_{h}^{i,T}$ is related to the index set of nodes with a strong connection to node *i*. In addition, coarse(*i, j, K _{h}*) is an appropriate cut-off (coarsening) function, e.g.,

$$\text{coarse}(i,j,{K}_{h})\u2254\alpha \cdot \underset{j,j\ne i}{\mathrm{max}}\{\mid {K}_{h}^{\left[ij\right]}\mid \},$$

(20)

with *α* [0, 1] (see, e.g., [22, §4.6.1]).

With those definitions a splitting into coarse and fine grid nodes can be achieved. For our application, a modified splitting algorithm is used [22, §4.6] as shown in Algorithm 2. Therein, the function

$$i\leftarrow \text{Pick}({\omega}_{h}\backslash ({\omega}_{C}\cup {\omega}_{F}\left)\right)$$

returns a node *i* for which the number $\mid {S}_{h}^{i,T}\mid +\mid {S}_{h}^{i,T}\cap {\omega}_{F}\mid $ is maximal. Note that tissue conductivity inhomogeneity and anisotropy are taken into account within the coarsening algorithm.

To achieve prolongation, the operator *P _{h,H}* :

$${P}_{h,H}^{\left[ij\right]}=\{\begin{array}{cc}1\hfill & i=j\in {\omega}_{C},\hfill \\ 1\u2215\mid {S}_{h}^{i,T}\cap {\omega}_{C}\mid \hfill & i\in {\omega}_{F},j\in {S}_{h}^{i,T}\cap {\omega}_{C},\hfill \\ 0\hfill & \text{else}.\hfill \end{array}\phantom{\}}$$

(21)

After the proper definition of the prolongation and coarse grid operators, it is possible to create in a recursive way a matrix hierarchy and an associated multigrid cycle, shown in Algorithm 3. Therein, the variable COARSEGRID denotes the level at which a direct solver is applied. For an m-*V* (*ν _{F}*,

The numerical examinations of the theory presented above were carried out in a four-layer sphere model with anisotropic skull compartment whose parameterization is shown in Table 1. For the choice of these parameters, we closely followed [12,15]. Forward solutions were computed for dipoles of 1 nAm amplitude located on the *y* axis at depths of 0% to 98.7% (in 1 mm steps) of the brain compartment (78 mm radius) using both radial (directed away from the center of the model) and tangential (directed parallel to the scalp surface) dipole orientations. *Eccentricity* is defined here as the percent ratio of the distance between the source location and the model midpoint divided by the radius of the inner sphere (78 mm). The most eccentric source position considered was thus only 1 mm below the CSF compartment. To achieve error measures which were independent of the specific choice of the sensor configuration, we distributed 748 electrodes in a regular fashion over the outer sphere surface. All simulations ran on a Linux-PC with an Intel Pentium 4 processor (3.2GHz) using the SimBio software environment [31] in which the algebraic multigrid solver package PEBBLES was embedded [9,38,42,10].

De Munck and Peters [5] derived series expansion formulas for a mathematical dipole in a multi-layer sphere model, denoted here as the *analytical solution*. The model consists of *S* shells with radii *r _{S}* <

$${\varphi}_{\mathrm{ana}}({x}_{0},{x}_{e})=\frac{1}{4\pi}\langle \mathbf{M},{S}_{0}\frac{{x}_{e}}{{r}_{e}}+({S}_{1}-\mathrm{cos}{\omega}_{0e}{S}_{0})\frac{{x}_{0}}{{r}_{0}}\rangle $$

with *ω*_{0e} the angular distance between source and electrode, and with

$${S}_{0}=\frac{{F}_{0}}{{r}_{0}}\frac{\Lambda}{{(1-2\Lambda \phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\omega}_{0e}+{\Lambda}^{2})}^{3\u22152}}+\frac{1}{{r}_{0}}\sum _{n=1}^{\infty}\left\{\right(2n+1\left){R}_{n}\right({r}_{0},{r}_{e})-{F}_{0}{\Lambda}^{n}\}{P}_{n}^{\prime}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\omega}_{0e}\right)$$

(22)

and

$${S}_{1}={F}_{1}\frac{\Lambda \phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\omega}_{0e}-{\Lambda}^{2}}{{(1-2\Lambda \phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\omega}_{0e}+{\Lambda}^{2})}^{3\u22152}}+\sum _{n=1}^{\infty}\left\{\right(2n+1\left){R}_{n}^{\prime}\right({r}_{0},{r}_{e})-{F}_{1}n{\Lambda}^{n}\}{P}_{n}\left(\mathrm{cos}\phantom{\rule{thinmathspace}{0ex}}{\omega}_{0e}\right).$$

(23)

The coefficients *R _{n}* and their derivatives, ${R}_{n}^{\prime}$, are computed analytically and the derivative of the Legendre polynomials,

$${t}_{k}\u2215{t}_{0}\le \upsilon ,\phantom{\rule{1em}{0ex}}{t}_{k}\u2254(2k+1){R}_{k}^{\prime}-{F}_{1}k{\Lambda}^{k}.$$

(24)

In the following simulations, a value of 10^{−6} was chosen for *ν* in (24). Using the asymptotic expansion, no more than 30 terms were needed for the series computation at each electrode.

The FE meshes of the four-layer sphere model were generated by the software TetGen [28] which used a *Constrained Delaunay Tetrahedralization* (CDT) approach [30,29]. This meshing procedure starts with the preparation of a suitable boundary discretization of the model in which for each of the layers and for a given triangle **edge** length, nodes are distributed in a regular fashion and connected through triangles. This yields a valid triangular surface mesh for each of the layers. Meshes of different layers are not intersecting each other. The CDT approach is then used to construct a tetrahedralization conforming to the surface meshes. It first builds a Delaunay tetrahedralization starting with the vertices of the surface meshes. The CDT then uses a local degeneracy removal algorithm combining vertex perturbation and vertex insertion to construct a new set of vertices which includes the input set of surface vertices. In a last step, a fast facet recovery algorithm is used to construct the CDT [30,29].

This approach is combined with two further constraints to the size and shape of the tetrahedra. The first constraint is important for the generation of quality tetrahedra. If *R* denotes the radius of the unique circumsphere of a tetrahedron and *L* its shortest edge length, the so-called *radius-edge ratio* of the tetrahedron can be defined as

$$\text{radius}-\text{edge}-\text{ratio}=R\u2215L.$$

(25)

The radius-edge ratio can distinguish almost all badly-shaped tetrahedra except one type of tetrahedra, so-called *slivers*. A sliver is a very flat tetrahedron which has no small edges, but can have arbitrarily large dihedral angles (close to *π*). For this reason, an additional mesh smoothing and optimization step is required to remove the slivers and improve the overall mesh quality.

A second constraint can be used to restrict the volume of the generated tetrahedra in a certain compartment. We follow the formula for regular tetrahedra:

$$\text{volume}=\sqrt{2}\u221512\cdot {\text{edge}}^{3}$$

(26)

Table 2 shows the number of **nodes** and **elements** of the six tetrahedra models used for the solver run-time comparison and accuracy tests. **factor** indicates the ratio of the number of nodes of the most highly resolved to both other models within each group. Additionally, the table contains the chosen **radius-edge-ratio** (see Equation (25)), the average **edge** length of the four triangular surface meshes, the corresponding **volume** constraints (see Equation (26)) for the tetrahedra and the compartments where the volume constraint is not applied. The most highly resolved meshes **tet503K** and **tet508K** of both groups had approximately the same resolution, while the others were chosen to have a factor of 4 coarser resolution with regard to the number of nodes. The meshes of group 1 concentrated the nodes in the outer three compartments because no volume constraint was applied for the inner brain compartment, while the nodes in the meshes of group 2 were distributed in a regular way throughout all four compartments. The meshes of group 1 were thus preferentially beneficial to the full subtraction approach, since the entries of the volume integral in Equation (10) are zero ((*σ*(*x*)–*σ*_{0}) = 0 for all *x* in the brain compartment) so that a coarse resolution can be expected to have no impact on the overall numerical accuracy, but will reduce the computational cost. In contrast, the meshes of group 2 were beneficial to both direct potential approaches. Figure 1 shows samples from the six tetrahedra models that were generated using the parametrizations from Table 2.

We compared numerical solutions with analytical solutions using three common error criteria [16,3,34,15,26]. The *relative (Euclidean) error* (RE) is defined as

$$\mathrm{RE}\u2254\frac{{\Vert {\underset{\u2012}{\varphi}}_{\mathrm{num}}-{\underset{\u2012}{\varphi}}_{\mathrm{ana}}\Vert}_{2}}{{\Vert {\underset{\u2012}{\varphi}}_{\mathrm{ana}}\Vert}_{2}},$$

where ${\underset{\u2012}{\varphi}}_{\mathrm{ana}}$, ${\underset{\u2012}{\varphi}}_{\mathrm{num}}\in {\mathbb{R}}^{\mathrm{m}}$ denote the analytical and the numerical solution vectors, respectively, at the m = 748 measurement electrodes. We furthermore defined

$$\mathrm{RE}(\%)\u2254100\cdot \mathrm{RE},\phantom{\rule{1em}{0ex}}\mathrm{max}\mathrm{RE}(\%)\u2254\underset{j}{\mathrm{max}}\left(\mathrm{RE}{(\%)}_{j}\right)$$

(27)

where *j* is the source eccentricity. In order to better distinguish between the topography (driven primarily by changes in dipole location and orientation) and the magnitude error (indicating changes in source strength), Meijs et al. [16] introduced the *relative difference measure* (RDM) and the *magnification factor* (MAG), respectively. For the RDM, we can show that

$$\mathrm{RDM}\u2254{\Vert \frac{1}{{\Vert {\underset{\u2012}{\varphi}}_{\mathrm{ana}}\Vert}_{2}}{\underset{\u2012}{\varphi}}_{\mathrm{ana}}-\frac{1}{{\Vert {\underset{\u2012}{\varphi}}_{\mathrm{num}}\Vert}_{2}}{\underset{\u2012}{\varphi}}_{\mathrm{num}}\Vert}_{2}=\sqrt{2\left(1-\mathrm{cos}\angle ({\underset{\u2012}{\varphi}}_{\mathrm{ana}},{\underset{\u2012}{\varphi}}_{\mathrm{num}})\right)}.$$

(28)

It therefore holds that 0 ≤ RDM ≤ 2, so that we can furthermore define

$$\mathrm{RDM}(\%)\u2254100\cdot \mathrm{RDM}\u22152,\phantom{\rule{1em}{0ex}}\mathrm{max}\mathrm{RDM}(\%)\u2254\underset{j}{\mathrm{max}}\left(\mathrm{RDM}{(\%)}_{j}\right).$$

(29)

The MAG is defined as

$$\mathrm{MAG}\u2254{\Vert {\underset{\u2012}{\varphi}}_{\mathrm{num}}\Vert}_{2}\u2215{\Vert {\underset{\u2012}{\varphi}}_{\mathrm{ana}}\Vert}_{2}$$

so that error minimum is at MAG = 1 and we therefore defined

$$\mathrm{MAG}(\%)=\mid 1-\mathrm{MAG}\mid \cdot 100,\phantom{\rule{1em}{0ex}}\mathrm{max}\mathrm{MAG}(\%)\u2254\underset{j}{\mathrm{max}}\left(\mathrm{MAG}{(\%)}_{j}\right).$$

(30)

With maxRE(%)^{k} we denote the maximal relative error in percent over all source eccentricities for an accuracy level of ACCURACY = 10^{−k}. The so-called *plateau-entry* for an iterative solver is then defined as the first *k* at which the condition

$$\mid \mathrm{max}\mathrm{RE}{(\%)}^{k}-\mathrm{max}\mathrm{RE}{(\%)}^{k+1}\mid \u2215\mathrm{max}\mathrm{RE}{(\%)}^{k+1}<0.05$$

(31)

is true.

The parameters of the Venant approach were chosen as proposed in [4]: The maximal dipole order (*n*_{0} in [4, Equation (22)]) and the scaling reference length (*α*_{ref} in [4, Equation (23)]) were set to *n*_{0} = 2 and *α*_{ref} = 20.0 mm, respectively. Since the chosen mesh size was a large factor smaller than the reference length, the second order term ${\left(\right(\Delta {\stackrel{\u2012}{x}}_{ki}^{r})}^{2}$ in [4, Equation (23)]) was small and the model focused on fulfilling the dipole moments of the zeros and first order. The exponent of the source weighting matrix (*n _{s}* in [4, Equation (25)]) was fixed to

The initial solution guess for all solvers was a zero potential vector. For the IC(0), ς = 0 was chosen for (14). For the AMG-CG, the 1-*V* (1, 1)-cycle AMG-preconditioner was used with *α* = 0.01 for (20). The factorization in Algorithm 3 was carried out whenever the size of the coarsest grid (COARSEGRID) in the preconditioner-setup was below 1000 and the coarse system was solved using a Cholesky-factorization. The setup times for the preconditioners were neglected in all calculations of computational cost because this step must be performed only once per head model. The evaluation with regard to relative solver accuracy in Algorithm 1 was limited to the discrete set of *accuracy levels *ACCURACY = 10^{−k} with k {0, … , 9}.

In a first study, we compared the numerical accuracy of the full subtraction approach (Section 2.1.1) with the two direct methods: Venant (Section 2.1.3) and partial integration (Section 2.1.2). Figure 2 shows the *RE*(%) for the different source eccentricities for the two finest models **tet503K** of group 1 (left) and **tet508K** of group 2 (right) (see Figure 1 and Table 2) with regard to the full subtraction (top row), the Venant (middle row) and the partial integration approach (bottom row). The results were computed with the AMG-CG and the necessary ACCURACY in Algorithm 1 for the plateau-entry (31) is indicated for both source orientation scenarios. In Figure 3, the maximal RE, RDM and MAG errors over all source eccentricities at the AMG-CG plateau-entry (31) are shown for all tetrahedra models, both source orientation scenarios and the three dipole modeling approaches.

Figure 2 clearly presents the advantages of the full subtraction approach whose error curves are smooth, while Venant and partial integration show an oscillating behavior. With RDM and MAG errors below 1% over all source eccentricities and for both orientation scenarios (see Figure 3), the full subtraction approach performs best for all source eccentricities for model **tet503K** (its mesh resolution was sufficiently high and the FE nodes were concentrated in the compartments CSF, skull and skin), where both direct approaches showed oscillations with a relatively high magnitude. As the results for model **tet508K** show, the oscillation magnitudes for the direct approaches could be strongly reduced by means of distributing the FE nodes in a regular way over all four compartments, hence decreasing the mesh size in the brain compartment. Nevertheless, even for model **tet508K**, the full subtraction approach was the most accurate method for nearly all source eccentricities. It was only outperformed by partial integration for the source which was only 1 mm below the CSF compartment. As both Figures Figures22 and and33 show, the partial integration approach performed well if the mesh was sufficiently fine in the brain compartment. The oscillation magnitudes of the Venant approach were generally even slightly smaller than for the partial integration approach, with only one exception (the result for the radial source 1 mm below the CSF compartment, shown in the middle row of Figure 2). The main reason for the outlier was that for the source 1 mm below the CSF, monopoles were positioned in the CSF compartment, which had a strong effect on the MAG for the radially oriented source.

Figure 4 shows the numerical error maxRE(%) versus the PCG solver ACCURACY from Algorithm 1 for the discrete set of accuracy levels from 10^{0} to 10^{−9}. Results for the high-resolution model **tet503K** of group 1 are shown in the left and from the high-resolution model **tet508K** of group 2 in the right column for the AMG-CG (top row), the IC(0)-CG (middle row) and the Jacobi-CG (bottom row).

maxRE(%) versus PCG solver ACCURACY (see Algorithm 1 and Section 3.5) for models **tet503K** of group 1 (left column) and **tet508K** of group 2 (right column) for the AMG-CG (top row), the IC(0)-CG (middle row) and the Jacobi-CG (bottom row). Source orientations **...**

The PCG accuracy measures the error in the solution vector of the FE linear equation system (8) (correction potential), (11) and (13) (total potential). For the full subtraction approach, maxRE(%) was thus not equal to 100 for ACCURACY = 10^{0} because ${\underset{\u2012}{\varphi}}_{\mathrm{num}}$ is equal to the analytically computed singularity potential ${\underset{\u2012}{\Phi}}_{0}$ from Equation (5). Because the PCG accuracy is measured in the ${K}_{h}{C}_{h}^{-1}{K}_{h}\text{-energy}$ norm, the plateau-entry (31) differs for different preconditioners *C _{h}*. As shown in Figure 4 for the high-resolution models and as collected in Table 3 for all six tetrahedra models, the maximally needed

In a last study, we compared solver wall-clock time versus numerical accuracy for the three different CG preconditioners AMG, IC(0) and Jacobi.

The average setup time for the AMG- and the IC(0)-preconditioner is given in Table 4. In the following, the time for the setup of the preconditioner was not included, because this step was carried out only once per head model.

In Figure 5, the solver time is shown versus the maxRE(%) for different levels of PCG accuracy for models **tet503K** and **tet33K** of group 1. The largest examined PCG ACCURACY level 10^{−k} is indicated in the figure. Please note that this level does not necessarily correspond to the plateau-entry level. In most cases results are presented up to one level higher.

Solver time versus maxRE(%) for models **tet503K** and **tet33K** of group 1 for tangentially and radially oriented sources for the potential approaches full subtraction (left), Venant (middle), and partial integration (right). Results are presented for the three **...**

For all tetrahedra models of groups 1 and 2, average solver times and iteration counts over all source eccentricities, source orientations and potential approaches for a plateau-entry (31) are collected in Table 5. Both Figure 5 and Table 5 clearly show the superiority of the AMG preconditioner. In all cases, even for the low-resolution grids **tet33K** and **tet32K**, the AMG-CG was the fastest solver, followed by the IC(0)-CG and the Jacobi-CG. The main result of Table 5 is the so-called *gain factor*, which is defined here as the result (solver time or iteration count) for the Jacobi-CG divided by the result for the AMG-CG. The gain factors clearly showed that the higher the mesh-resolution, i.e., the higher the condition number of the corresponding FE stiffness matrix, the larger the difference in performance between AMG-CG, IC(0)-CG, and Jacobi-CG. An increasing mesh-resolution led to a strong increase in the number of iterations of IC(0)-CG (factor of 3.2 between **tet503K** and **tet33K** and 4.1 between **tet508K** and **tet32K**) and Jacobi-CG (factor of 3.0 between **tet503K** and **tet33K** and 3.6 between **tet508K** and **tet32K**), while the number of AMG-CG iterations was only slightly increasing (factor of 1.9 between **tet503K** and **tet33K** and 1.8 between **tet508K** and **tet32K**). This clearly shows the stronger *h*-dependence of the IC(0) and Jacobi preconditioners.

The goals of this technical study of finite element (FE) based solution techniques for the electroencephalographic forward problem were twofold. The first aim was to compare three efficient iterative FE solver techniques under realistic conditions that still allowed quasi-analytical solutions. The second aim was to evaluate three different numerical formulations of the current dipole, which is the bioelectric source most commonly used to represent neural electrical activity. A major motivation of such studies is the special need to achieve high accuracy and efficiency with FE based approaches for this problem. The many advantages of this approach are often hindered by the unacceptable computation costs of implementing it so that improved efficiency will provide substantial progress to the field.

When using the ${K}_{h}{C}_{h}^{-1}{K}_{h}\text{-energy}$ norm stopping criterion for the PCG algorithm applied on meshes with up to 500K nodes, a relative solver accuracy of 10^{−6} for AMG-CG, 10^{−7} for IC(0)-CG and 10^{−8} for Jacobi-CG was necessary and sufficient to fall below the discretization error. The AMG-CG achieved an order of magnitude higher computational speed than the CG with the standard preconditioners with an increasing gain factor with decreasing mesh size. The increasing gain factor shows that the convergence rate of the Jacobi- and IC(0)-preconditioning methods are much more dependent on the mesh size *h* than the AMG-CG, as discussed in detail in Section 2.2.3. However, while for the geometric multigrid, an *h*-independent convergence rate and an *h*-independent condition number can be proven ([11, Lemma 10.7.1,Theorem 10.7.15]), the AMG-CG was not optimal in our application with a slight *h*-dependence shown by a slightly increasing iteration count with increasing mesh resolution. Such a result had to be expected because the source analysis stiffness matrix was not an M-matrix and the prolongation operator of the presented AMG-CG was tuned for speed and not for an optimal behavior with regard to the iteration count. A discrete harmonic extension as proposed in [22] improved the interpolation properties, but the application of this prolongation operator is more expensive, which decreased the overall run-time performance in our application.

We generated two groups of Constrained Delaunay tetrahedralization (CDT) FE meshes [30,29], tuned for the specific needs of the different potential approaches. In group 1, for the full subtraction approach [7], FE nodes were concentrated in the CSF, skull and skin, while the brain compartment was meshed as coarsely as possible. Group 2 was tuned for the needs of both direct potential approaches [45,4,1,37], which profit more from a regular distribution of FE nodes over all four compartments and especially a higher resolution at the source positions.

With regard to the numerical error, in the tuned FE meshes with about 500K nodes we achieved high accuracies—in the range of a few percent maximal relative error (maxRE)—over all source eccentricities for both the full subtraction and the two direct potential approaches. With a maximal relative difference measure (maxRDM) and a maximal magnification factor (maxMAG) of less than 1% over all source eccentricities for sources up to 1 mm below the CSF compartment (model **tet503K**, maximal examined eccentricity of 98.7%), the full subtraction approach performed consistently better than both direct approaches. We found that the instability of the full subtraction approach with regard to the RE for high source eccentricities was mainly a magnitude instability (MAG), not a topographic one (RDM). Our results clearly illustrate the advantages of the full subtraction approach as long as the homogeneity condition is sufficiently fulfilled, i.e., as long as the distance of the source to the next conductivity inhomogeneity is large enough or the resolution of the FE mesh at the nearest conductivity inhomogeneity to the source is fine enough. A theoretical reasoning for this finding is given in [39]. While error curves oscillated for both direct approaches, they were smooth for the full subtraction approach as long as the homogeneity condition was sufficiently fulfilled. The smoothness can be explained by the fact that the singularity potential is sensitively modeling the details of the point dipole source singularity while the FE method is able to accurately model the smooth correction potential. The oscillating behavior of the direct approaches can easily be understood for the partial integration approach. Since linear basis functions were used, the right-hand side vector (12) and thus the computed solution for the total potential in (11) is constant over any finite element for all sources with moment **M**_{0} and any location **x**_{0} within the finite element. For partial integration with linear basis functions, the inverse source localization resolution is thus bounded by the size of the element, which should be fine as long as high-resolution FE meshes are used. The direct approach oscillations should therefore not necessarily be a problem with regard to the inverse problem. It has furthermore been shown using quasi-analytically computed reference potentials in a four-layer anisotropic sphere model and single dipole fit reconstructions based on the Venant FE approach in a regular tetrahedral FE mesh with 161K nodes, that localization errors were in the 1 millimeter range [40] and might thus most often be neglected when compared to errors due to, e.g., EEG data noise, inaccurate tissue segmentation and conductivity modeling.

Schimpf et al. [26] investigated different FE potential approaches in a four-layer sphere model with isotropic skull and sources up to 1 mm below the CSF compartment. In their report, a regular 1 mm cube model was used (thus a much higher FE resolution) and a maxRDM of 7% and a maxMAG of 25% was achieved with a subtraction approach, which performed best in their comparison. Awada et al. [1] implemented a two-dimensional subtraction approach and compared its numerical accuracy with a partial integration method in a two-dimensional multi-layer sphere model. A direct comparison with our results is therefore difficult, but the authors concluded that the subtraction method was more accurate than the partial integration direct approach. In a locally refined (around the source singularity) tetrahedral mesh with 12,500 nodes of a four-layer sphere model with anisotropic skull and first order FE basis functions in a subtraction approach, Bertrand et al. [3] reported a maxRDM of above 20% and a maxMAG up to 70% for a maximal eccentricity of 97.6%. Van den Broek [34] used a subtraction approach in a locally refined (around the source singularity) tetrahedral mesh with 3,073 nodes of a three-layer sphere model with anisotropic skull. For the maximal examined eccentricity of 94.2%, they reported a maxRDM of up to 50%.

However, the right-hand side (RHS) vector is expensive to compute and is densely populated (i.e., *N _{h}* non-zeros) for the full subtraction approach (10) and sparse with just some few (|NODESOFELE| for partial integration (12), and

The following limitations of the presented work are important and can be seen as an outlook for our future studies. Oscillations of the direct potential approaches might be avoided. For example, the use of higher order basis functions might cure the oscillation problem for the partial integration approach. For both examined direct methods, the oscillations might be avoided by means of precomputing a lead field (i.e., many forward solutions) for those source positions with minimal errors (e.g., only at the barycenters of brain finite elements for the partial integration approach and only at brain FE nodes for the Venant approach) and consequently using a lead field interpolation technique [46] for all other forward solutions. Lead field extrapolation [46] has to be examined as a concept to further decrease the numerical error for eccentric sources. Such an approach might be particularly fruitful for the subtraction method, where a violation of the homogeneity condition led to especially large numerical errors or instabilities as shown in Figure 2 for the most eccentric source in especially the model **tet508K** of group 2. In our study, we examined numerical errors for a single source at each discrete 1mm step eccentricity level through the inner compartment along the *y* axis. By means of the error oscillations and magnitudes for the direct methods and the specific sensitivity of the subtraction method for eccentric sources, the dependence of all three dipole modeling approaches on the specific FE mesh properties has been shown. A further important measure for establishing the validity of the presented results would be a more statistical error measure where, for each eccentricity level, not only the error of a single source on the *y* axis is evaluated, but an error statistic is evaluated over many sources with different positions in the volume conductor and/or sources with slightly varying (e.g., 0.1 mm) source positions around each eccentricity level.

The AMG-CG turned out to achieve an order of magnitude higher computational speed than Jacobi-CG or incomplete Cholesky-CG for the FEM based EEG forward and inverse problem. Our results corroborate the theoretical results that the higher the FE resolution, the greater the advantage of using MG preconditioning. The AMG-CG together with the fast transfer matrix approach now enable resolutions which seemed to be impracticable before. In the comparison of dipole modeling approaches, highest accuracies were achieved with the full subtraction approach in CDT meshes where nodes were concentrated in the compartments CSF, skull and skin. However, the accuracy of the subtraction approach is especially dependent on a sufficiently fulfilled homogeneity condition, so that the direct potential approaches partial integration and Venant, which also performed well in our comparison, might be even more appropriate if white matter conductivity anisotropy is taken into account.

This research was supported by the German Research Foundation (DFG), projects WO1425/1-1 and JU445/5-1, and by the Center for Integrative Biomedical Computing, NIH NCRR Project 2-P41-RR12553-07. The authors would like to thank S. Reitzinger, G. Haase and H. Si for providing the PEBBLES algebraic multigrid code and the Tetgen FE meshing tool, resp., and for their quick responses whenever needed and the anonymous reviewers for their helpful critics and comments that significantly improved our manuscript.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

[1] Awada K, Jackson D, Williams J, Wilton D, Baumann S, Papanicolaou A. Computational aspects of finite element modeling in EEG source localization. IEEE Trans Biomed. Eng. 1997;44(8):736–751. [PubMed]

[2] Baillet S, Mosher J, Leahy R. Electromagnetic brain mapping. IEEE Signal Processing Magazine. 2001 November;18(6):14–30.

[3] Bertrand O, Thévenet M, Perrin F. 3D finite element method in brain electrical activity studies. In: Nenonen J, Rajala H, Katila T, editors. Biomagnetic Localization and 3D Modelling. 1991. Report of the Dep. of Tech.Physics, Helsinki University of Technology.

[4] 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. Electroenc. Clin. Neurophysiol. 1997;102:267–278. [PubMed]

[5] de Munck J, Peters M. A fast method to compute the potential in the multi sphere model. IEEE Trans Biomed. Eng. 1993;40(11):1166–1174. [PubMed]

[6] de Munck J, van Dijk B, Spekreijse H. Mathematical dipoles are adequate to describe realistic generators of human brain activity. IEEE Trans Biomed. Eng. 1988;35(11):960–966. [PubMed]

[7] Drechsler F, Wolters C, Dierkes T, Si H, Grasedyck L. A highly accurate full subtraction approach for dipole modelling in EEG source analysis using the finite element method. Max-Planck-Institut für Mathematik in den Naturwissenschaften; 2007. Techn. Report 95. http://www.mis.mpg.de/publications/preprints/2007/prepr2007-95.html.

[8] Gencer N, Acar C. Sensitivity of EEG and MEG measurements to tissue conductivity. Phys.Med.Biol. 2004;49:701–717. [PubMed]

[9] Haase G, Kuhn M, Reitzinger S. Parallel AMG on distributed memory computers. SIAM J. Sci.Comp. 2002;24(2):410–427.

[10] Haase G, Reitzinger S. Cache issues of algebraic multigrid methods for linear systems with multiple right-hand sides. SIAM J. Sci.Comp. 2005;27(1):1–18.

[11] Hackbusch W. Iterative solution of large sparse systems of equations, Springer Verlag. Applied Mathematical Sciences. 1994;95

[12] Hallez H, Vanrumste B, Hese PV, D’Asseler Y, Lemahieu I, de Walle RV. A finite difference method with reciprocity used to incorporate anisotropy in electroencephalogram dipole source localization. Phys.Med.Biol. 2005;50:3787–3806. [PubMed]

[13] Jung M, Langer U. Applications of multilevel methods to practical problems. Surv.Math.Ind. 1991;1:217–257.

[14] Kickinger F. Algebraic multigrid for discrete elliptic second-order problems. In: Hackbusch W, editor. Multigrid Methods V. Proceedings of the 5th European Multigrid conference; Berlin. Springer; 1998.

[15] Marin G, Guerin C, Baillet S, Garnero L, Meunier G. Influence of skull anisotropy for the forward and inverse problem in EEG: simulation studies using the FEM on realistic head models. Human Brain Mapping. 1998;6:250–269. [PubMed]

[16] Meijs J, Weier O, Peters M, van Oosterom A. On the numerical accuracy of the boundary element method. IEEE Trans Biomed. Eng. 1989;36:1038–1049. [PubMed]

[17] Mohr M. Simulation of Bioelectric Fields: The Forward and Inverse Problem of Electroencephalographic Source Analysis. vol. Band 37 of Arbeitsberichte des Instituts für Informatik. Friedrich-Alexander-Universität Erlangen-Nürnberg; 2004. iSSN 1611-4205.

[18] Murakami S, Okada Y. Contributions of principal neocortical neurons to magnetoencephalography and electroencephalography signals. J.Physiol. 2006;575(3):925–936. [PubMed]

[19] Ollikainen J, Vauhkonen M, Karjalainen P, Kaipio J. Effects of local skull inhomogeneities on EEG source estimation. Med.Eng.Phys. 1999;21:143–154. [PubMed]

[20] Plonsey R, Heppner D. Considerations on quasi-stationarity in electro-physiological systems. Bull.Math.Biophys. 1967;29:657–664. [PubMed]

[21] Reitzinger S. Ph.D. thesis. Schriften der Johannes-Kepler-Universität Linz, Reihe C - Technik und Naturwissenschaften; 2001. Algebraic multigrid methods for large scale finite element equations. No.36.

[22] Ruge J, Stüben K. Algebraic multigrid (AMG) In: McCormick S, editor. Multigrid Methods. vol. 5 of Frontiers in Applied Mathematics. SIAM; Philadelphia: 1986. pp. 73–130.

[23] Rullmann M, Anwander A, Dannhauer M, Warfield S, Duffy F, Wolters C. EEG source analysis of epileptiform activity using a 1mm anisotropic hexahedra finite element head model. NeuroImage. 2009;44(2):399–410. http://dx.doi.org/10.1016/j.neuroimage.2008.09.009. [PMC free article] [PubMed]

[24] Saad Y. Iterative methods for sparse linear systems. PWS Publishing Company; 1996.

[25] Sarvas J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Phys.Med.Biol. 1987;32(1):11–22. [PubMed]

[26] Schimpf P, Ramon C, Haueisen J. Dipole models for the EEG and MEG. IEEE Trans Biomed. Eng. 2002;49(5):409–418. [PubMed]

[27] Schwarz H. Methode der finiten Elemente. B.G.Teubner Stuttgart; 1991.

[28] Si H. TetGen, a quality tetrahedral mesh generator and three-dimensional delaunay triangulator, v1.3, user’s manual. Weierstrass Institute for Applied Analysis and Stochastics; 2004. Tech. Rep. 9.

[29] Si H. Adaptive tetrahedral mesh generation by constrained Delaunay refinement. International Journal for Numerical Methods in Engineering. 2008;75(7):856–880.

[30] Si H, Gärtner K. Meshing piecewise linear complexes by constrained Delaunay tetrahedralizations. Proc. 14th International Meshing Roundtable; Sandia National Laboratories; 2005.

[31] SimBio -IPM and -NeuroFEM: Algorithms for inverse EEG and MEG source analysis. 2000. ( http://www.simbio.de) https://www.mrt.uni-jena.de/neurofem/index.php/Main_Page. since.

[32] Stüben K. A review of algebraic multigrid. Nov., 1999. Techn. Report 69, GMD, http://www.gmd.de.

[33] Tanzer O, Järvenpää S, Nenonen J, Somersalo E. Representation of bioelectric current sources using whitney elements in finite element method. Phys.Med.Biol. 2005;50:3023–3039. [PubMed]

[34] van den Broek S. Ph.D. thesis. Proefschrift Universiteit Twente Enschede; The Netherlands: 1997. Volume conduction effects in EEG and MEG.

[35] Waberski T, Buchner H, Lehnertz K, Hufnagel A, Fuchs M, Beckmann R, Rienäcker A. The properties of source localization of epileptiform activity using advanced headmodelling and source reconstruction. Brain Top. 1998;10(4):283–290. [PubMed]

[36] Weinstein D, Zhukov L, Johnson C. Lead-field bases for EEG source imaging. Annals of Biomed. Eng. 2000;28(9):1059–66. [PubMed]

[37] Weinstein D, Zhukov L, Johnson C. Lead-field bases for electroencephalography source imaging. Annals of Biomed.Eng. 2000;28(9):1059–1066. [PubMed]

[38] Wolters C. Influence of Tissue Conductivity Inhomogeneity and Anisotropy on EEG/MEG based Source Localization in the Human Brain. Leipzig: 2003. (No. 39 in MPI Series in Cognitive Neuroscience, MPI of Cognitive Neuroscience Leipzig). iSBN 3-936816-11-5 (also: Univ., Diss., http://lips.informatik.uni-leipzig.de/pub/2003-33/en)

[39] Wolters C, Köstler H, Möller C, Härtlein J, Grasedyck L, Hackbusch W. Numerical mathematics of the subtraction method for the modeling of a current dipole in EEG source reconstruction using finite element head models. SIAM J. on Scientific Computing. 2007;30(1):24–45.

[40] Wolters CH. Habilitation in Mathematics, Faculty of Mathematics and Natural Sciences. University of Münster; Germany: 2008. Finite element method based electro- and magnetoencephalography source analysis in the human brain.

[41] Wolters CH, Anwander A, Berti G, Hartmann U. Geometry-adapted hexahedral meshes improve accuracy of finite element method based EEG source analysis. IEEE Trans.Biomed.Eng. 2007;54(8):1446–1453. [PubMed]

[42] Wolters CH, Anwander A, Reitzinger S, Haase G, Halgren E, Ahlfors S, Hämäläinen M, Cohen D. Algebraic multigrid with multiple right-hand-side treatment for an efficient computation of EEG and MEG lead field bases. BIOMAG2004, Proc. of the 14th Int. Conf. on Biomagnetism; Boston, USA. Aug.8-12,, 2004; http://www.biomag2004.org.

[43] Wolters CH, Grasedyck L, Hackbusch W. Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem. Inverse Problems. 2004;20(4):1099–1116.

[44] Wolters CH, Köstler H, Möller C, Härdtlein J, Anwander A. Numerical approaches for dipole modeling in finite element method based source analysis. International Congress Series 1300; June 2007; pp. 189–192. iSBN-13:978-0-444-52885-8, http://dx.doi.org/10.1016/j.ics.2007.02.014.

[45] Yan Y, Nunez P, Hart R. Finite-element model of the human head: Scalp potentials due to dipole sources. Med.Biol.Eng.Comput. 1991;29:475–481. [PubMed]

[46] Yvert B, Crouzeix-Cheylus A, Pernier J. Fast realistic modeling in bioelectromagnetism using lead-field interpolation. Human Brain Mapping. 2001;14:48–63. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |