Elsevier

NeuroImage

Volume 46, Issue 4, 15 July 2009, Pages 1055-1065
NeuroImage

A full subtraction approach for finite element method based source analysis using constrained Delaunay tetrahedralisation

https://doi.org/10.1016/j.neuroimage.2009.02.024Get rights and content

Abstract

A mathematical dipole is widely used as a model for the primary current source in electroencephalography (EEG) source analysis. In the governing Poisson-type differential equation, the dipole leads to a singularity on the right-hand side, which has to be treated specifically. In this paper, we will present a full subtraction approach where the total potential is divided into a singularity and a correction potential. The singularity potential is due to a dipole in an infinite region of homogeneous conductivity. The correction potential is computed using the finite element (FE) method. Special care is taken in order to evaluate the right-hand side integral appropriately with the objective of achieving highest possible convergence order for linear basis functions. Our new approach allows the construction of transfer matrices for fast computation of the inverse problem for anisotropic volume conductors. A constrained Delaunay tetrahedralisation (CDT) approach is used for the generation of high-quality FE meshes. We validate the new approach in a four-layer sphere model with a highly conductive cerebrospinal fluid (CSF) and an anisotropic skull compartment. For radial and tangential sources with eccentricities up to 1 mm below the CSF compartment, we achieve a maximal relative error of 0.71% in a CDT-FE model with 360 k nodes which is not locally refined around the source singularity and therefore useful for arbitrary dipole locations. The combination of the full subtraction approach with the high quality CDT meshes leads to accuracies that, to the best of the author's knowledge, have not yet been presented before.

Introduction

Inverse methods are used to reconstruct current sources in the human brain by means of electroencephalography (EEG) or magnetoencephalography (MEG) measurements of, e.g., event related fields or epileptic seizures (Hämäläinen et al., 1993, Michel et al., 2004, Plummer et al., 2008, Rullmann et al., 2009). A critical component of the inverse neural source reconstruction is the solution of the forward problem, i.e., the simulation of the fields at the head surface for a known primary current source in the brain. Because of the availability of quasi-analytical EEG forward problem solution formulas (de Munck and Peters, 1993), the head volume conductor is still often represented by a multi-layer sphere model. However, this model is just a rough approximation to the reality, so that numerical approximation methods are more and more frequently used such as the boundary element method (BEM) (Hämäläinen et al., 1993, Kybic et al., 2005), the finite volume method (FVM) (Cook and Koles, 2006), the finite difference method (FDM) (Hallez et al., 2005) or the finite element method (FEM). This paper focuses on the FEM which allows high accuracy for the numerical solution of elliptic partial differential equations since it is specifically tailored to the corresponding variational formulation (Hackbusch, 1992, Braess, 2007, Dahmen and Reusken, 2008) and since it allows high flexibility in modelling the forward problem in geometrically complicated inhomogeneous and anisotropic head volume conductors (Bertrand et al., 1991, Awada et al., 1997, van den Broek, 1997, Marin et al., 1998, Schimpf et al., 2002, Zhang et al., 2006, Wolters et al., 2007a, Wolters et al., 2007b). Advantages of the FEM are that the underlying weak formulation of the boundary value problem and its Galerkin-discretisation allow on the one side an appropriate and accurate modelling of the physics and offer on the other side a mathematically clean treatment of cases with low regularity (Hackbusch, 1992, Braess, 2007, Dahmen and Reusken, 2008) like the one in EEG source analysis with the discontinuous tissue conductivities. As will be shown, the FEM discretisation can be adapted to the structure of the solution function by means of appropriate mesh refinements and coarsenings so that highest accuracy can be achieved with a reasonable number of nodes.

It was shown in Sarvas (1987), de Munck et al. (1988), Hämäläinen et al. (1993) and Murakami and Okada (2006) that the mathematical dipole is an adequate model to represent the primary current which is caused by a synchronous activity of tens of thousands of densely packed apical dendrites of large pyramidal cells oriented in parallel in the human cortex. The dipole model is thus considered to be the “atomic” structure of the primary current density distribution that has to be reconstructed within the inverse problem. Hence, one of the key questions for all 3D forward modelling techniques is the appropriate modelling of the potential singularity introduced into the differential equation by means of the mathematical dipole.

Direct potential approaches (Yan et al., 1991, Buchner et al., 1997) approximate the dipole moment through optimally distributed monopolar sources and sinks on neighbouring FE nodes of the source location. This approach leads to finite distances between the poles that seem reasonable and it performs well in validation studies (Buchner et al., 1997, Wolters et al., 2007b). However, a disadvantage of direct approaches is the absence of a well-understood mathematical theory. Furthermore, in recent comparison studies of different direct methods with the subtraction approach (Awada et al., 1997, Schimpf et al., 2002), it is concluded that the overall best accuracy is achieved using the latter method.

A subtraction approach for the modelling of a mathematical dipole in FE-based source analysis was widely suggested (Bertrand et al., 1991, Awada et al., 1997, van den Broek, 1997, Marin et al., 1998, Schimpf et al., 2002, Wolters et al., 2007a). All proposed approaches have in common that the total potential is divided into an analytically known singularity potential and a singularity-free correction potential which can then be approximated numerically using an FE approach. The subtraction approach is applicable for arbitrary geometries and the FE discretisation allows for a treatment of geometric details and tissue inhomogeneities and anisotropies. In Wolters et al. (2007a), a theoretical insight was given into the subtraction approach. A proof was given for existence and uniqueness of a weak solution in the function space of zero-mean potential functions and convergence properties of the FE-approach to the correction potential were stated. A projected subtraction method was proposed where the singularity potential was projected in the FE space (Wolters et al., 2007a). This approach was shown to perform well in a three-compartment (skin, skull, brain) sphere model provided that the so-called source eccentricity was limited to 95%. The eccentricity is generally defined as the percent ratio of the distance between the source location and the model midpoint divided by the radius of the inner sphere. When considering a three-shell model, 95% eccentricity seems reasonable because the dipoles that are located in the cortex will have an eccentricity even lower than 92% as reported in Marin et al. (1998).

However, the three-compartment model of the head, which is still most often used in source analysis (see, e.g., Hämäläinen et al., 1993, Michel et al., 2004, Kybic et al., 2005, Wolters et al., 2007a), ignores the cerebrospinal fluid (CSF) compartment between the cortex and the skull. (Baumann et al., 1997) showed for a group of 7 human subjects (neurosurgical patients, three males, four females, ranging in age from 4.5 months to 70 years) that the CSF conductivity at body-temperature can accurately be measured to be 1.79 S/m with a maximal standard deviation of less than 1.4% between subjects for the large frequency range between 10 Hz and 10,000 Hz. Additionally, the CSF compartment was shown to have a significant influence on EEG source analysis (Ramon et al., 2004, Wendel et al., 2008, Rullmann et al., 2009). In four-compartment models, this layer is taken into account, but source eccentricity then has to be determined with regard to the inner CSF surface, i.e., the most eccentric sources are only 1 or 2 mm apart from the next conductivity discontinuity. Therefore, eccentricities of more than 98% have to be examined. It is well-known (and in Wolters et al., 2007a, a theoretical reasoning was given for this fact), that with increasing eccentricity, the numerical accuracy decreases. This is not only the case for the subtraction approach (Bertrand et al., 1991, van den Broek, 1997, Marin et al., 1998, Schimpf et al., 2002, Wolters et al., 2007a), but also for the direct approaches in FE modelling (Yan et al., 1991, Buchner et al., 1997, Schimpf et al., 2002, Wolters et al., 2007b). Examining modelling accuracy at high source eccentricities is thus important and constitutes a major numerical challenge. In this paper, a four-compartment sphere model is used for validation of the proposed FE approach, the CSF compartment is modelled with the isotropic conductivity value of 1.79 S/m (Baumann et al., 1997) and sources are considered with eccentricities up to 98.7% (1 mm below the CSF). Furthermore, as Marin et al. (1998) and Wolters et al. (2007a) have shown, modelling skull conductivity anisotropy is a further numerical challenge and generally leads to higher numerical errors than modelling the skull as an isotropic compartment. Following Marin et al. (1998), we will model the skull compartment with a conductivity anisotropy of 1:10.

The numerical accuracy of a subtraction approach is determined by means of the following two important criteria: i) FE approach for the numerical computation of the correction potential and ii) FE mesh generation procedure. Until now, in 4 layer sphere model validation studies with realistic source eccentricities up to 98%, none of the presented FE approaches were shown to achieve satisfying numerical accuracy in combination with being practically feasible in an inverse source analysis scenario. The projected subtraction approach presented in Wolters et al. (2007a) only performed well up to 95% eccentricity. Schimpf et al. (2002) investigated an FE subtraction approach in a regular 1 mm hexahedral model of a four layer sphere volume conductor with isotropic skull and sources up to 1 mm below the CSF compartment. For the most eccentric source, a topography error of 7% and a magnitude error of 25% were reported. It is well-known also from other studies that the abrupt transitions and right angles of regular hexahedral models at material interfaces negatively influence especially the magnitude error (Wolters et al., 2007b). Even if a geometry-adaptation of the hexahedra can alleviate this problem (Wolters et al., 2007b), a well-chosen tetrahedralisation will still better represent smooth tissue boundaries and can furthermore be adapted to the specific structure of the solution function. However, until now, only ordinary Delaunay tetrahedralisation as defined in Si and Gärtner (2005) and Si (2008) was proposed for FE-based source analysis (Bertrand et al., 1991, van den Broek, 1997, Buchner et al., 1997, Marin et al., 1998, Wolters et al., 2007a). In Bertrand et al. (1991), van den Broek (1997), and Marin et al. (1998), coarse ordinary tetrahedral meshes were considered yielding unacceptably large numerical errors already at eccentricities above 90%. In Bertrand et al. (1991) and van den Broek (1997), local mesh refinement around the source was used to achieve better results. However, with regard to the inverse problem, the setup of source-location dependent locally refined meshes is difficult to implement and time-consuming to compute and thus might not be practicable for an inverse source analysis.

In this paper, we propose a so-called full subtraction FE approach which appropriately evaluates the right-hand side integral for the correction potential with the objective of achieving highest possible convergence order for linear basis functions. A constrained Delaunay tetrahedralisation (CDT) approach will be presented and used for the generation of high-quality FE meshes. Even if we will also propose to adapt the tetrahedralisation by means of an appropriate refinement and coarsening to the correction potential solution function, our new approach does not need local mesh refinement around the source. As we will show, it therefore allows the construction of transfer matrices for fast computation of the inverse problem. We validate the new approach in a four-layer sphere model with highly-conductive CSF and anisotropic skull compartment and sources up to 1 mm below the CSF compartment. We compare the accuracy of our new method with the projected subtraction approach from Wolters et al. (2007a) and the literature. It will be shown that the combination of the full subtraction approach with the CDT-FE meshes leads to very low numerical errors for the computed EEG potentials. Note that since the magnetoencephalography (MEG) forward problem is also based on the computed electric potential (see, e.g., Hämäläinen et al., 1993, Wolters et al., 2004), our method is directly applicable to MEG source analysis, too.

Section snippets

The continuous forward problem

The mathematical model for the numerical simulation of electric and magnetic fields in the human head is based on the quasistatic approximation of Maxwell's equations (Plonsey and Heppner, 1967, Sarvas, 1987, Hämäläinen et al., 1993). Following Sarvas (1987), de Munck et al. (1988), Hämäläinen et al. (1993) and Murakami and Okada (2006), the primary current density function is described byJp=Jy,Jy(x)=divM(y)δ(xy),yYΩ,M(y)R3with M(y) being the current dipolar moment at position y, ΩR3 the

Evaluation with regard to right-hand side quadrature order

In the first study, we compared the numerical accuracy of the presented full subtraction approach for quadrature formulas with different integration order for the right-hand side (Eqs. (12) and (13)). The goal of this study was to verify that second order integration formulas are necessary and sufficient as stated in the section “Finite element approach”. Fig. 2 shows the relative errors between the numerical and the quasi-analytical solutions for tangential (left column) and radial sources

Discussion and conclusion

We presented theory and numerical experiments of a full subtraction approach to model a mathematical dipole in finite element (FE) method based electroencephalography (EEG) source reconstruction using constrained Delaunay tetrahedral (CDT) meshes. Since the magnetoencephalography (MEG) forward problem is also based on the computed electric potential (see, e.g., (Hämäläinen et al., 1993, Wolters et al., 2004)), our method is directly applicable to MEG source analysis. We embedded the approach

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (WO1425/1-1, GR3179/1-1, JU 445/5-1). The authors would like to thank J.C. de Munck for providing the software for the quasi-analytical solution in multilayer sphere models and for his quick responses whenever needed.

References (46)

  • de MunckJ. et al.

    A fast method to compute the potential in the multi sphere model

    IEEE Trans. Biomed. Eng.

    (1993)
  • de MunckJ. et al.

    Mathematical dipoles are adequate to describe realistic generators of human brain activity

    IEEE Trans. Biomed. Eng.

    (1988)
  • GencerN. et al.

    Sensitivity of EEG and MEG measurements to tissue conductivity

    Phys. Med. Biol.

    (2004)
  • HaaseG. et al.

    Parallel AMG on distributed memory computers

    SIAM J. Sci. Comput.

    (2002)
  • HackbuschW.

    Elliptic Differential Equations. Theory and Numerical Treatment

    (1992)
  • HackbuschW.
  • HallezH. et al.

    A finite difference method with reciprocity used to incorporate anisotropy in electroencephalogram dipole source localization

    Phys. Med. Biol.

    (2005)
  • HämäläinenM. et al.

    Magnetoencephalography: theory, instrumentation, and applications to noninvasive studies of the working human brain

    Rev. Mod. Phys.

    (1993)
  • KybicJ. et al.

    A common formalism for the integral formulations of the forward eeg problem

    IEEE Trans. Med. Imag.

    (2005)
  • MarinG. et al.

    Influence of skull anisotropy for the forward and inverse problem in EEG: simulation studies using the FEM on realistic head models

    Hum. Brain Mapp.

    (1998)
  • MeijsJ. et al.

    On the numerical accuracy of the boundary element method

    IEEE Trans. Biomed. Eng.

    (1989)
  • MurakamiS. et al.

    Contributions of principal neocortical neurons to magnetoencephalography and electroencephalography signals

    J. Physiol.

    (2006)
  • PlonseyR. et al.

    Considerations on quasi-stationarity in electro-physiological systems

    Bull. Math. Biophys.

    (1967)
  • Cited by (0)

    View full text