736 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 8, AUGUST 1997 Computational Aspects of Finite Element Modeling in EEG Source Localization Kassem A. Awada, Member, IEEE, David R. Jackson,* Senior Member, IEEE, Jeffery T. Williams, Member, IEEE, Donald R. Wilton, Fellow, IEEE, Stephen B. Baumann, Member, IEEE, and Andrew C. Papanicolaou Abstract— A comparison is made of two different implemen- tations of the finite element method (FEM) for calculating the potential due to dipole sources in electroencephalography (EEG). In one formulation (the direct method) the total potential is the unknown that is solved for and the dipole source is directly incorporated into the model. In the second formulation (the subtraction method) the unknown is the difference between the total potential and the potential due to the same dipole in an infinite region of homogeneous conductivity, corresponding to the region where the dipole is located. Both methods have the same FEM system matrix. However, the subtraction method requires an additional calculation of flux integrations along the edges of the elements in the computation of the right-hand side (RHS) vector. It is shown that the subtraction method is usually more accurate in the forward modeling, provided the flux integrations are computed accurately. Errors in calculating the flux integrations may result in large errors in the forward solution due to the ill-conditioned nature of the FEM system matrix caused by the Neumann boundary condition. To minimize the errors, closed-form expressions for the flux integrations are used for both linear and quadratic triangular elements. It is also found that FEM forward modeling errors may cause false extrema in the least-square objective function obtained from the boundary potential, near boundaries between media of differing conductivity. Multiple initial guesses help eliminate the possibility of the solution getting trapped in these false extrema. Index Terms—electroencephalography, finite element method, inverse problems, imaging, numerical methods. I. INTRODUCTION M OST of the source localization techniques used in electroencephalography (EEG) are based on models assuming one or more dipole sources. The forward solution for the potential is first calculated from the dipole model, and the dipole parameters (location and orientation) are then found by minimizing the least-square error between the calculated and measured boundary potential. Since errors in the modeling of potential will translate into errors in source localization, it Manuscript received December 27, 1995; revised March 7, 1997. This work was supported in part by the State of Texas Advanced Technology Program. Asterisk indicates corresponding author. K. A. Awada and S. B. Baumann are with the Department of Neurological Surgery, University of Pittsburgh, Pittsburgh, PA 15213 USA. *D. R. Jackson is with the Department of Electrical and Computer En- gineering, University of Houston, Houston, TX 77204-4793 USA (e-amil: djackson@uh.edu). J. T. Williams and D R. Wilton are with the Department of Electrical and Computer Engineering, University of Houston, Houston, TX 77204-4793 USA. A. C. Papanicolaou is with the Department of Neurosurgery, The University of Texas—Houston Health Science Center, Houston, TX 77030 USA. Publisher Item Identifier S 0018-9294(97)05353-6. is important to have an accurate calculation of the forward solution. Realistic models of the geometry can improve the accuracy of forward solutions. Head shapes can be obtained from magnetic resonance images (MRI’s) by extracting surface boundaries for the major tissues, such as scalp, skull, cere- brospinal fluid (CSF) and brain. The boundary element method (BEM) assumes that tissues are homogeneous and isotropic between the boundaries. Several studies have shown that the BEM can improve the accuracy of EEG and magnetoen- cephalographic (MEG) source localization, when compared to the commonly used homogeneous sphere or concentric-sphere models of the head ([1]–[6]). However, most tissues in the brain are not homogeneous nor isotropic. To account for tissue heterogeneity more tissues than the four common ones (scalp, skull, CSF, and brain), as well as tissue subtypes such as hard or soft bone, need to be extracted from MRI scans. Some recent studies using the finite element method (FEM) have used as many as 12 tissues, plus air in the sinuses, and conductivities were assigned using values in the literature ([7] and [8]). Although anisotropy was not accounted for in these studies, finite element models have the capability of assigning a conductivity tensor to each element and are, therefore, potentially more accurate. Some tissues, such as white matter in the brain and muscle in the head, are known to be anisotropic [9], but a full accounting of anisotropy for all the tissues in the human head has yet to be performed. That should not deter the development of finite element models of the head, itself a daunting task, that can later be refined with more accurate and complete values for tissue con- ductivity and anisotropy. There are a number of fundamental questions that must be answered regarding the implementation of finite element models for source localization in the brain. For example, what is the best way to implement the forward solution? How coarse can the mesh be and under what circumstances does it need to be locally refined? Are there circumstances under which some tissues can be ignored? How important is anisotropy and how much does it affect solution accuracy? Can parallelized algorithms speed computation to the point of having practical turn-around times for solutions? How can multiple or diffuse sources be handled efficiently? Three-dimensional (3-D) FEM studies have recently been performed that attempt to address some of these issues ([7], [8], and [10]). In this paper we address one of these questions by presenting a direct and detailed comparison between forward solutions obtained by two distinct implementations of the FEM that 0018–9294/97$10.00  1997 IEEE