PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
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.006
PMCID: PMC2791331
NIHMSID: NIHMS120338

Accuracy and run-time comparison for different potential approaches and iterative solvers in finite element method based EEG source analysis

Abstract

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.

Keywords: electroencephalography, source reconstruction, finite element method, dipole singularity, full subtraction potential approach, Venant potential approach, partial integration potential approach, preconditioned conjugate gradient method, algebraic multigrid, incomplete Cholesky, Jacobi, constrained Delaunay tetrahedralization, anisotropic four-layer sphere model

1. Introduction

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

2. Theory

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 jp is governed by the Poisson equation with homogeneous Neumann boundary conditions on the head surface Γ = [partial differential]Ω [20,25], which we can express as

(σΦ)=jp=JpinΩ,σΦ,n=0onΓ,
(1)

with n the unit surface normal, and assuming a reference electrode with given potential, i.e., Φ(xref) = 0. The primary currents are modeled by a mathematical dipole at position x0R3 with moment M0R3 [25,6,18],

Jp=jp(x)=(M0δ(xx0)),
(2)

where δ is the Dirac delta distribution.

2.1. Finite element modeling techniques for the potential singularity

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.

2.1.1. Full subtraction approach

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

Φ=Φ0+Φcorr,
(3)

where the singularity potential, Φ0, is defined as the solution for a dipole in an unbounded homogeneous conductor with constant conductivity σ0. σ0R3×3 is the conductivity at the source position, which is assumed to be constant in a non-empty subdomain Ω0 around x0, in the following called the homogeneity condition. The solution of Poisson’s equation under these conditions for the singularity potential

(σ0Φ0)=jp
(4)

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

Φ0(x)=14πdetσ0M0,(σ0)1(xx0)(σ0)1(xx0),(xx0)32.
(5)

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

(σΦcorr)=((σ0σ)Φ0)inΩ,
(6)

with inhomogeneous Neumann boundary conditions at the surface:

σΦcorr,n=σΦ0,nonΓ.
(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 [var phi]i. When projecting the correction potential into the FE space, i.e., Φcorr(x)Φcorr,h(x)=j=1Nhφj(x)ucorr,h[j] (Nh being the number of FE nodes), and applying variational and FE techniques to (6) and (7), we finally arrive at a linear system [7]

Khucorr,h=jcorr,h,
(8)

with the stiffness matrix

Kh[i,j]=Ωσφj,φidx,
(9)

for KhRNh×Nh, and the right-hand side vector jcorr,hRNh with entries

jcorr,h[i]=Ω(σ0σ)Φ0,φi(x)dxΩφi(x)n(x),σ0Φ0(x)dx.
(10)

We then seek for the coefficient vector ucorr,h=(ucorr,h[1],,ucorr,h[Nh])RNh 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].

2.1.2. The partial integration direct approach

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

Ω(σΦ)φidx=Ωjpφidx.

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

ΩσΦ,φidx+ΓσΦ,nφidΓ=Ωjp,φidx+Γjp,nφidΓ.

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

ΩσΦ,φidx=Ωjp,φidx=(2)M0,φi(x0).

Setting Φ(x)Φh(x)=j=1Nhφj(x)uh[j], leads to the linear system

Khuh=jPI,h,
(11)

with the same stiffness matrix as in (9) and the right-hand side vector jPI,hRNh with entries

jPI,h[i]={M0,φi(x0)ifiNODESOFELE(x0),0otherwise.}
(12)

The function NODESOFELE(x0) determines the set of nodes of the element which contains the dipole at position x0. Note that while the right-hand side vector (10) is fully populated, jPI,h has only |NODESOFELE| non-zero entries. Here, |·| denotes the number of elements in the set NODESOFELE. For the linear basis functions [var phi]i considered here, the right-hand side (12) and thus the computed solution for the total potential in (11) will be constant for all x0 within a finite element.

2.1.3. The Venant direct approach

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

Khuh=jVenant,h
(13)

with the same stiffness matrix as in (9). The right-hand side vector jVenant,hRhN 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].

2.2. FE solver methods

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 Kh is equal to

κ2(Kh)=λmaxλmin

with λmax the largest and λmin the smallest eigenvalues, respectively, of Kh [11, §2.10]. The condition number behaves asymptotically as O(h2) and condition numbers of more than 107 have been computed for FE problems in EEG source analysis [38]. Large condition numbers are the reason for slow convergence of common iterative solvers [11,24] and any effective solution approach has to minimize the effects of this poor conditioning.

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 jh of the linear equation system [11, §3.4]. The goal of a preconditioner, ChRNh×Nh, is the reduction of κ2(Ch1Kh) for the preconditioned equation system Ch1Khuh=Ch1jh. Further requirements are that it is cheap with regard to arithmetic and memory costs to solve linear systems Chwh = rh with rh the residual and wh the residual for the preconditioned system.

Algorithm 1
PCG: (Kh, uh, jh, Ch, ACCURACY) → (uh)

Theorem 2.1 (Error estimate for PCG method) Let Kh and Ch be positive definite. If uh denotes the exact solution of the equation system, then the kth iterate of the PCG method uhk fulfills the following energy norm estimate

uhkuhKhck=21+c2kuh0uhKh,cκ2(Ch1Kh)1κ2(Ch1Kh)+1.

Proof: Hackbusch [11, Theorem 9.4.14].

As indicated in Algorithm 1, the PCG method is stopped after the kth iteration if the relative error, i.e., ehk=uhkuh in the controllable KhCh1Kh-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 Kh with a vector and the preconditioning step Chwh = rh. We used the compact row storage format for sparse matrices so that both the storage of the stiffness matrix and the matrix-vector multiplication are of complexity O(Nh) (see [27, Chapter 4.6.4],[24,10]).

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

2.2.1. Jacobi preconditioning or scaling

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

ChDh2,DhDIAG(Kh[11],,Kh[NhNh]).

When splitting the Jacobi-preconditioner between left and right (row and column scaling), one has to solve K~hvh=Dh1jh with K~h=Dh1KhDhtr and uh=Dhtrvh. Row and column scaling preserves symmetry, so that the scaled matrix K~h is again SPD with unit diagonal entries. The scaling may lead to a first substantial condition improvement, because diagonal entries in Kh 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 Kh be SPD and ChDh2 the Jacobi-preconditioner. Assume that each row of Kh does not contain more than d nonzero entries. Then, for all diagonal matrices D~h1, it is

κ2(Ch1Kh)dκ2(D~h1Kh),

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

Proof: Hackbusch [11, Theorem 8.3.3].

2.2.2. Incomplete Cholesky preconditioning

The SPD stiffness matrix Kh can be decomposed into a left triangular matrix Lh and its transpose using the Cholesky-decomposition, Kh=LhLhtr [27, pp.209f]. Nevertheless, because of a large fill-in, ChLhLhtr would not be appropriate as a preconditioner. The Incomplete Cholesky (IC) preconditioner without fill-in, IC(0), is defined as ChL0L0tr where L0 is the Cholesky-decomposition of the scaled stiffness matrix K~h which is restricted to the same non-zero-pattern as the lower triangular part of K~h. For incomplete factorizations, the preconditioning operation Chwh = rh in Algorithm 1 is solved by a forward-back sweep, which can efficiently use the compact row storage format for the sparse matrix L0 as described in detail in [27, Chapter 4.6.4]. The existence of IC(0) is not necessarily guaranteed for general SPD matrices. Therefore, a reduction of non-diagonal stiffness matrix entries has to be carried out in certain applications before IC(0) computation is possible [27, p.266]. If the scaled stiffness matrix is decomposed by means of K~h=Eh+Idh+Ehtr, with E~hRNh×Nh its strict lower triangular part, the reduction can be formulated as

K˘h=Idh+11+ς(Eh+Ehtr).
(14)

For sufficiently large ςR0+, the existence of IC(0) is guaranteed, but with increasing ς, the preconditioning effect decreases. Note that for certain special cases, a condition improvement to O(h1) can be proven as, e.g., when using a modified ILUω-preconditioning with ω = −1 (in the symmetric case, the ILU0 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]).

2.2.3. Algebraic multigrid preconditioning

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

κ2(Ch1Kh)11ρm,
(15)

with Ch the preconditioner resulting from m steps of the GMG method. As shown in [13,11,32], a robust method which provides a small convergence rate for a wide class of real-life problems is given by exploiting the MG-method as a preconditioner for the CG method. With MG(m)-CG, we denote the MG-preconditioned CG method with m the number of MG iterations for the CG preconditioning step. The GMG(m)-CG can improve the convergence rate to ρ/4, if ρ is assumed to be small, as shown in [11, §10.8.3].

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 Kh (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 ith row of Kh is considered as being related to a grid point in ωh (the index set of nodes), and an off-diagonal entry is related to an edge in an FE grid. A description of AMG is now given for a symmetric two grid method, where h is related to the fine grid and H to the coarse grid. Each AMG algorithm consists of the following components:

  1. Coarsening: define the splitting ωh = ωC [union operator] ωF of ωh into sets of coarse and fine grid nodes ωC and ωF, respectively.
  2. Transfer operators: prolongation Ph,H:RNHRNh and its adjoint as the restriction
    RH,hPh,Htr.
    (16)
  3. Definition of the coarse matrix by Galerkin’s method, i.e.,
    KHRH,hKhPh,H.
    (17)
    Because of (b), KHRNH×NH is again SPD.
  4. Appropriate smoother for the considered problem class: In order to achieve a symmetric method, e.g., a forward Gauss-Seidel method for pre-smoothing and the adjoint, a backward Gauss-Seidel method for post-smoothing ([11, §4.8.3,§10.7.1,2],[22, §4.4]).

Coarsening

The coarsening process has the task of reducing the number of nodes such that NH = |ωC| < Nh = |ωh|. The grid points ωh can be split into two disjoint subsets ωC (coarse grid nodes) and ωF (fine grid nodes), i.e., ωh = ωC [union operator]ωF and ωC [intersection operator]ωF = [empty] such that there are (almost) no direct connections between any two coarse grid nodes and such that the resulting number of coarse grid nodes is as large as possible [32, p.12]. Instead of considering all connections between nodes as being of the same rank, the following sets are introduced

Nehi={jKh[ij]ζKh[i,i],ij},
(18)
Shi={jNehiKh[ij]>coarse(i,j,Kh)},Shi,T={jNehiiShj)},
(19)

where Nehi is the index set of neighbors (a pre-selection is carried out by the threshold-parameter ςR0+), Shi denotes the index set of nodes with a strong connection from node i and Shi,T is related to the index set of nodes with a strong connection to node i. In addition, coarse(i, j, Kh) is an appropriate cut-off (coarsening) function, e.g.,

coarse(i,j,Kh)αmaxj,ji{Kh[ij]},
(20)

with α [set membership] [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

iPick(ωh\(ωCωF))

returns a node i for which the number Shi,T+Shi,TωF is maximal. Note that tissue conductivity inhomogeneity and anisotropy are taken into account within the coarsening algorithm.

Algorithm 2
COARSE:({Shi,T},ωh)(ωC,ωF)

Prolongation

To achieve prolongation, the operator Ph,H : VH [mapsto] Vh has to be defined correctly. The form that turned out to be the most efficient for the presented application was proposed in [14] and is given by

Ph,H[ij]={1i=jωC,1Shi,TωCiωF,jShi,TωC,0else.}
(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, νB)-cycle AMG preconditioned CG method, the operation SOLVE Chwh = rh in Algorithm 1 is realized by m calls of MG(Kh,wh, rh, νF, νB). Specific performance optimizations for the AMG preconditioner used in this study are described in detail in [10, Section 2].

Algorithm 3
V-cycle MG : (Kh, uh, jh, νF, νB) → (uh)

3. Methods

3.1. Validation platform

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

Table 1
Parameterization of the anisotropic four-layer sphere model.

3.2. Analytical solution in an anisotropic multilayer sphere model

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 rS < rS−1 < … < r1 and constant radial, σrad(r)=σjradR+, and constant tangential conductivity, σtang(r)=σjtangR+, within each layer rj+1 < r < rj. It is assumed that the source at position x0 with radial coordinate r0R is in a more interior layer than the measurement electrode at position xeR3 with radial coordinate re=r1R. The spherical harmonics expansion for the mathematical dipole (2) is expressed in terms of the gradient of the monopole potential to the source point. Using an asymptotic approximation and an addition-subtraction method to speed up the series convergence yields

ϕana(x0,xe)=14πM,S0xere+(S1cosω0eS0)x0r0

with ω0e the angular distance between source and electrode, and with

S0=F0r0Λ(12Λcosω0e+Λ2)32+1r0n=1{(2n+1)Rn(r0,re)F0Λn}Pn(cosω0e)
(22)

and

S1=F1Λcosω0eΛ2(12Λcosω0e+Λ2)32+n=1{(2n+1)Rn(r0,re)F1nΛn}Pn(cosω0e).
(23)

The coefficients Rn and their derivatives, Rn, are computed analytically and the derivative of the Legendre polynomials, Pn, are determined by means of a recursion formula. We refer to [5] for the derivation of the above series of differences and for the definition of F0, F1 and [turned v]. Here, it is only important that the latter terms are independent of n and that they can be computed from the given radii and conductivities of layers between source and electrode and of the radial coordinate of the source. The computations of the series (22) and (23) are stopped after the k-th term if the following criterion is fulfilled

tkt0υ,tk(2k+1)RkF1kΛ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.

3.3. Tetrahedral mesh generation

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

radiusedgeratio=RL.
(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:

volume=212edge3
(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.

Fig. 1Fig. 1Fig. 1Fig. 1Fig. 1Fig. 1
Cross-sections of the six tetrahedral meshes of the four compartment sphere model. The corresponding parametrizations of the models are shown in Table 2. Visualization was done using the software TetView [28].
Table 2
The six tetrahedra models used for the solver time comparison and accuracy tests. The table shows the number of nodes and elements of each mesh and factor indicates the ratio of the number of nodes of the most highly resolved to both other models within ...

3.4. Error criteria

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

REϕnumϕana2ϕana2,

where ϕana, ϕnumRm denote the analytical and the numerical solution vectors, respectively, at the m = 748 measurement electrodes. We furthermore defined

RE(%)100RE,maxRE(%)maxj(RE(%)j)
(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

RDM1ϕana2ϕana1ϕnum2ϕnum2=2(1cos(ϕana,ϕnum)).
(28)

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

RDM(%)100RDM2,maxRDM(%)maxj(RDM(%)j).
(29)

The MAG is defined as

MAGϕnum2ϕana2

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

MAG(%)=1MAG100,maxMAG(%)maxj(MAG(%)j).
(30)

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

maxRE(%)kmaxRE(%)k+1maxRE(%)k+1<0.05
(31)

is true.

3.5. FEM and solver parameter settings

The parameters of the Venant approach were chosen as proposed in [4]: The maximal dipole order (n0 in [4, Equation (22)]) and the scaling reference length (αref in [4, Equation (23)]) were set to n0 = 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 ((Δxkir)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 (ns in [4, Equation (25)]) was fixed to ns = 1 and the regularization parameter (λD in [4, Equation (25)]) was chosen as λD = 10−6. These settings effect a spatial concentration of the monopole loads in the dipole axis around the dipole location.

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 = 10k with k [set membership] {0, … , 9}.

4. Results

4.1. Numerical error versus potential approach

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.

Fig. 2Fig. 2Fig. 2Fig. 2Fig. 2Fig. 2
RE(%) versus source eccentricity for the two most highly resolved models tet503K of group 1 (left) and tet508K of group 2 (right) using the full subtraction (top row), the Venant (middle row) and the partial integration (bottom row) potential approaches. ...
Fig. 3Fig. 3Fig. 3Fig. 3Fig. 3Fig. 3
maxRE(%), maxRDM(%) and maxMAG(%) accuracies for the full subtraction, the Venant and the partial integration approach for all six tetrahedra models (see Figure 1 and Table 2) and both tangential (left) and radial (right) source orientation scenarios ...

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.

4.2. Numerical error versus PCG accuracy

Figure 4 shows the numerical error maxRE(%) versus the PCG solver ACCURACY from Algorithm 1 for the discrete set of accuracy levels from 100 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).

Fig. 4Fig. 4Fig. 4Fig. 4Fig. 4Fig. 4
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 = 100 because ϕnum is equal to the analytically computed singularity potential Φ0 from Equation (5). Because the PCG accuracy is measured in the KhCh1Kh-energy norm, the plateau-entry (31) differs for different preconditioners Ch. As shown in Figure 4 for the high-resolution models and as collected in Table 3 for all six tetrahedra models, the maximally needed k (for a PCG accuracy of ACCURACY = 10k) decreased when the preconditioning quality increased (except for the radial source orientation in model tet503K, see Fig. 4). Furthermore, as Table 3 shows, a higher PCG accuracy was needed for the plateau-entry when the mesh resolution increased.

Table 3
Maximally needed k [set membership] {0, … , 9} for a PCG ACCURACY = 10k for the plateau-entry (31) over all three potential approaches.

4.3. Numerical error versus solver time

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.

Table 4
Average setup time (sec.) for the AMG- and the IC(0)-preconditioner.

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 10k 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.

Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5Fig. 5
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.

Table 5
Average solver time (sec.) and iteration count (iter) over all source eccentricities, source orientations and potential approaches for plateau-entry (31). For all tetrahedra models of groups 1 and 2, results are presented for the three different CG preconditioners ...

5. Discussion

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 KhCh1Kh-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 M0 and any location x0 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., Nh non-zeros) for the full subtraction approach (10) and sparse with just some few (|NODESOFELE| for partial integration (12), and C for Venant (13)) non-zero vector entries for the direct approaches, which has implications for the computational effort when using the fast FE transfer matrix approach for EEG and MEG [43] (additionally, see [37,8,12]), which limits the total number of FE linear equation systems to be solved for any inverse method to the number of sensors m. After solving m FE linear equation systems to compute the transfer matrix, each forward problem can be solved by a single multiplication of the RHS vector with the transfer matrix [43], resulting in a computational effort of 2 * m * P operations with P = Nh for the full subtraction, P = |NODESOFELE| for partial integration and P = C for the Venant approach. Note that the transfer matrix approach can not be used if the mesh is adapted according to varying source positions within the inverse problem. We therefore attempted to avoid local mesh refinement techniques as used in [3,34].

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.

6. Conclusion

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.

Acknowledgment

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.

Footnotes

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.

References

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