Realistic meshless conductor model for EEG inverse problems

International Journal of Bioelectromagnetism

http://www.ijbem.org

Vol. 10, No. 3, pp. 176 - 189, 2008

Realistic Meshless Conductor Model for EEG Inverse Problems Mauricio Pohl-Alfaro, Oscar Yáñez-Suárez, Juan R. Jimenéz-Alaniz, Verónica Medina-Bañuelos Neuroimaging Laboratory,Universidad Autónoma Metropolitana UAM-Iztapalapa, México D.F., México Correspondence: Mauricio Pohl, Laboratorio de Neuroimagenología, Departamento de Ingeniería Eléctrica, Universidad Autonóma Metropolitana- Iztapalapa, México D.F., México. E-mail: [email protected], phone +52 55 58044903

Abstract. In this work, a meshless (free Galerkin) method is used to define a realistic head conductor model to solve the EEG inverse problem. The aim is to find sources responsible for scalp potentials with precision and low computation time. In order to assess the proper method’s performance, it was evaluated with simulated EEG signals, obtained from different source positions and constant electrode number. The meshless method integrates anatomical information obtained from segmented anatomical MRI; its performance was compared with that obtained considering an infinite homogenous conductor model. Head geometry and secondary current sources, contributing to the potential recorded at electrode sites, are taken into account in the meshless method, thus obtaining a realistic head conductor volume model; this results in source localization errors near 13 ± 5 mm on 200 cases of simulated EEG records. Keywords: EEG inverse problem; meshless method; sLORETA; mean shift MRI segmentation;

1. Introduction Neural current sources produce external magnetic fields and electric potentials on scalp surface; these can be measured using magnetoencephalography (MEG) and electroencephalography (EEG), respectively (Mosher et al., 1999; Spinelli et al., 2000; Eskola, 1999). Head current fields producing these MEG and EEG signals can be separated into two components: a primary current term, representing the impressed neural and passive cellular currents, and the secondary or volume currents that are a result of the macroscopic electric field (Mosher et al., 1999; Tripp, 1983). Primary currents are the sources of interest in E/MEG, because they represent neuronal activity areas associated with cognitive, motor and sensorial processes of the brain (Srinivasan, 1999; Teubner et al., 2005; Neilson et al., 2005; Wu et al., 2003). Localization of the cortical regions responsible for this activity is of importance in the correct diagnosis and treatment of several brain pathologies; for instance, surgical treatment of epilepsy requires the exact localization of critical activity foci, during the preoperative planning phase, in order to avoid affecting eloquent areas (Mirkovic et al., 2003). The clinical standard used to obtain a larger precision in source localization is electrocorticography (ECoG), but due to the fact that it is a highly invasive, risky procedure, its use becomes limited (Zhang and van Drongelen, 2003). Thus, electric potentials recorded at scalp, provide a non-invasive means to monitor spatiotemporal evolution of brain neuronal activity (Ding et al., 2005, Darvas and Pantazis, 2004). In this case, EEG inverse solutions are proposed to estimate electrical activity sources, from those signals recorded at scalp; however, the main drawback is that standard 10-20 EEG does not contain enough information about spatial activity distribution (Tissari and Rahola, 2003; Zou et al., 2004; PascualMarqui and Michel, 1994) and therefore, from these data, a unique solution to the inverse problem cannot be found (Pascual-Marqui and Esslen, 2002). Several clinical techniques have been used to identify sources through an inverse procedure; functional neuroimages based on hemodynamics (PET, SPECT, fMRI) provide a good localization of

176

generators and do not require an a priori conductor model; techniques based on electromagnetic source localization (EEG, MEG) require a brain volume in which inverse solutions are solved and current sources estimated over a finite set of locations. Several brain conductor models can be found in the literature, going from concentric spheres with homogeneous conductivity or different sphere combinations (de Munck and Peters, 1993) to realistic individual models (He et al., 1987, Hämäläinen and Sarvas, 1989, Hedou, 2004; Ellenrieder et al., 2005; van `t Ent et al., 2001) , that take into account the head geometry and where the solution space can be restricted either to the cortex or to the whole volume; several comparisons have been carried out between these approaches to locate brain sources (Mosher et al. 1999; Zhai and Yao, 2004; Wu et al., 2003; Akalm et al., 2004; Cuffin, 1995; Rao et al., 2000). In contrast to the spherical model, several steps are required to build up a realistic head representation; first, anatomical information (MRI or CT) must be obtained from the patient; these images must be segmented in order to extract the surface and/or different brain structures where the solutions must be estimated. In this procedure, a robust segmentation method becomes necessary to adequately define gray matter as the solution space. Afterwards, this solution space must be geometrically described, either using a mesh or through a parametric description; for this purpose, four methods have been commonly employed: a) Finite differences b) Boundary Element Method (BEM), based on a reformulation of the problem as a surface integral equation c) Finite Element Method (FEM) where volume is described as tetrahedral elements, with the possibility to incorporate anisotropic characteristics and d) Projective methods that include the meshless or free-Galerkin method where electromagnetic fields are described as a combination of weighted basis functions (Belytschko et al., 1997; Xuan et al., 2004; Zhang et al., 2005). In both, FEM and BEM, node connections become a burdening, computationally expensive task; and also the methods need a special symmetrical collocation method for the nodes. For the meshless method, these drawbacks are overcome, since no connection between nodes is required. Recent research developed on meshless methods applied to the ECG direct problem, conclude that these methods are comparable in efficiency, complexity and accuracy to FEM (Li et al., 2007). The next step in the definition of a realistic head model needs the exact localization of electrode positions, so they can be correctly registered and labeled in the anatomical images, Finally, several numerical methods can be used to solve the direct problem directly into the real model, that is, finding the electrode potentials given a known source. This article introduces the lead field matrix needed to find the inverse solution of Poisson’s equation, considering a real, anatomically based conductor volume defined with a meshless method. This inverse formulation stems from the contribution of Ellenrieder (Ellenrieder et al., 2005, Ellenrieder et al., 2006) that propose a direct problem EEG solution, also considering a meshless model. In the next section, the theoretical background corresponding to the direct problem is presented, followed by the proposed inverse solution formulation; implementation and validation procedures are also described, including image segmentation, EEG signal simulation and inverse solution strategy (sLORETA). In section 3 results of the method’s performance are presented, to finally discuss and conclude in the last section.

2. Material and Methods 2.1. Theoretical formulation of Direct and Inverse Meshless Methods. Direct Meshless Method. Electric potential Φ(r ) at any location r within the solution space can be described as a superposition of two terms: Φ(r ) = Φ F (r ) + Φ N (r )

(1)

where the first term in Eq. 1 corresponds to primary current sources and the second one represents secondary or superficial sources. The associated direct problem consists on solving Poisson’s equations with Neumann boundary conditions to compute these potentials given a source (for r at the scalp, the problem is thus known as the direct EEG problem). It can be formulated as

177

∇ 2 Φ F (r )=∇ 2 Φ(r )=∇ ⋅(qδ (r − p )) ⎧⎪∇ 2 Φ N (r )= 0; (2) ⎨ ⎪⎩∇Φ N (r ) ⋅n δ Ω G = − ∇Φ F (r ) ⋅n δ Ω G where, q and p are magnitude and position of the current source confined in a space or body ΩG, and boundary δΩG that represents the subject’s head; the component ΦF(r) corresponds to electric potential at point r produced by the source, in an infinite homogeneous medium with electrical conductivity σ (Eq. 3):

Φ F (r ) =

1 (r − p ) ⋅q 4πσ r −p 3

(3)

The function ΦF(r), as well as its gradient (projected over the surface normal n) can be analytically computed and therefore the only term that must be numerically calculated is ΦN(r). In the estimation of the electric potential registered on the surface of a 3D conducting body of arbitrary shape and constant conductivity, generated by a source distribution inside it, the size and electromagnetic features of the body allow a quasi-static approximation of Maxwell equations (Hämäläinen and Sarvas, 1989). When formulating the boundary integral equation (BIE), nodes are placed in m locations ri , i=1..m, both inside and over the surface of the body, associating to each one an equation that describes the potential, thus obtaining a linear set of equations that represents the system; the solution is a vector of potentials in each one of the nodes. These equations relate the potential ΦN(ri) in each node ri, with nearby nodes rj , j = 1..n, which lie inside a region of influence (ROI) of node ri; each region is a subset of the total volume ( Ω L ⊂ Ω G ) and the union of all of them i corresponds to the total volume of the body ( ∪ΩL = ΩG ). i In the numerical method proposed by (Ellenrieder et al., 2005), the electrical potential is taken as a local linear combination (over an spherical ROI) of a set of weighted functions, such that, for any r ∈ ΩG , the electrical potential due to secondary currents is estimated as:

ˆ (r ) = pT (r )a Φ N (r ) ≈ Φ N

(4)

where p(r) is a vector of m approximation basis functions. In a simple case these basis functions can be assumed as order 2 polynomials, defined as (m = 10):

) [

(

p(r ) = p rx , ry , rz = 1, rx , ry , rz , rx ry , rx rz , ry rz , rx2 , ry2 , rz2

]

T

(5)

The vector a contains the linear combination coefficients, that are chosen in such a way that they minimize the estimation error in the Laplace equation solution in a spatial region ΩG. At the same time, they should satisfy the border conditions of the surface node set. Due to the local character of the approximation, the coefficient vector a varies in space; in order to find these coefficients it becomes necessary to use moving least squares (Eq. 6), so as to reduce the error between the potential ~ (r )) and the potential previously estimated by Eq. 4. This is coefficients assumed for each node ( Φ i N accomplished by minimizing the following Jacobian function: B(ri ) =

∑ w (r )(Φ n

~

j

i

N (ri ) − p

T

(r j )a

)

2

(6)

j =1

where n is the number of nodes within the spherical ROI belonging to node ri. In our case, the sphere radius is R and wj(r) are Gaussian weighting functions of neighboring node j, defined by: ⎧ − ( r − r j λ j )2 − (R λ j )2 −e ⎪e r − rj ≤ R ⎪ − (R λ j )2 (7) w j (r ) = ⎨ 1− e ⎪ r − rj ≥ R ⎪⎩0

178

where λj is a positive attenuation constant and R is the ROI radius. To estimate a unique solution to the moving least squares minimization problem, the number of nodes in the ROI of node r has to be greater than m. From the Jacobian minimization in Eq.6 we obtain the following solution to the approximation of the potential in node r: ~ ˆ (r ) = Ψ T (r )Φ Φ N N =

n

∑Ψ

~

(8)

j (r )Φ N j

j =1

where Ψj is the shape function of node j and is given by the following expression, which shows its dependence on the basis (P = [p(rj)]) and weighting functions (W(r) = [wj(r)]) vectors:

(

ΨT (r ) = pT (r) PT W(r )P

)

−1

PT W(r )

(9)

The approximation of the potential derivatives from Eq.8 turns out to be: ˆ (r ) ⋅n(r ) = ∇Φ N

n

∑∇Ψ

j (r ) ⋅n

(r ) Φ~ N j

(10)

j =1

It should be noted that Ψj(ri) ≠ δij and therefore the method does not interpolate over the electrical potential values of the nodes, as BEM or FEM do. The following linear system of equations is formulated by defining an equation for each node of the total volume and considering that the potential laplacian due to secondary sources is equal to zero and that the current flow normal to the surface can be analytically determined:

~ H ΦN = F

(11)

where vector F is composed by the potential gradients fk due to primary sources at the surface nodes rk which is defined by: f k = ∇ Φ N (rk ) ⋅n(rk ) = − ∇ Φ F (rk ) ⋅n(rk )

(12)

After solving Eq. 11, the estimated potential coefficients are given by:

~ Φ N = H −1F

(13)

Thus, the electric potential due to secondary sources for each node is defined as:

~ Φ N = ΨΦ N

(14)

Inverse Meshless Method. As mentioned in the previous section, the scalp potential can be described as the addition of two factors: the first is obtained by supposing that the current source is immersed in an infinite homogeneous medium (measured at scalp level) and the second is due to the secondary current density. The total potential at the EEG recording electrodes is given by an expression analogous to Eq.1:

Φ e = Φ eF + Φ eN

(15)

where Φe = [Φ(e1), Φ(e2),….Φ(eN)]T, is the vector of scalp potentials recorded at N electrodes placed at locations ei. The Poisson equation relating the current densities inside the volume with each electrode’s potential is given by: ∇ ⋅ (σ ∇Φ e ) = −∇ ⋅ J

(16)

179

where J = [J1,J2,…,JM]T, is the current density vector for M nodes inside region ΩG. The potential Φe due to the current density inside the volume, can be expressed as a matrix product, of the lead field matrix K, that holds the volume’s geometrical and electrical information, multiplied by the vector of current densities inside this volume. And therefore:

Φe = K J

(17)

Using Eq. 17 into Eq. 15 gives:

K J = ΦeF + ΦeN

(18)

The relationship between current source and potential ΦeF, can be found analytically. Let K ∞ be the lead field matrix in an infinite homogeneous medium; the following expression can then be defined:

Φ eF = K ∞ J

(19)

where the elements of K ∞ are given by:

(

)

k ∞ eβ = k ∞ re , rβ =

1 4πσ

re − rβ

(20)

3

re − rβ

re are electrode coordinates, where e =1…N rβ are voxel coordinates, with β = 1…M Substituting Eq.19 into Eq. 18, we obtain:

K J = K ∞ J + Φ eN

(21)

The term ΦeN corresponds to electric potential due to secondary sources at the electrodes. By incorporating Eq. 13 and Eq. 14 into Eq. 21, we get:

K J = K ∞ J + ΨH −1F (22) The vector F is given by the gradient’s potential projected over the normal surface, denoted by: F = ∇Φ eN ⋅ n

(23)

Because of Neumann boundary conditions, we know that the gradient potential due to a superficial charge on one side must be balanced with the volumetric charge inside the body. Thus, for those electrodes placed at the surface, we can write:

∇Φ eN ⋅ n ∂Ω = − ∇Φ eF ⋅ n ∂Ω L

(24)

L

By inserting Eq. 24 into Eq. 22, we obtain: K J = K ∞ J + ΨH −1 (− ∇Φ eF ⋅ n )

(25)

and from Eq. 19 we have: K J = K ∞ J − ΨH −1K ∞ J

(26)

Therefore, the lead field matrix K required to solve the inverse problem, taking into account a realistic head model generated with a meshless method, is given by:

180

(

K = K M = K ∞ I − ΨH −1

)

(27)

2.2 Implementation a) Inverse problem solution using standardized Low Resolution Electromagnetic Tomography (sLORETA) The estimates of the distributed current sources Jˆ can be found by solving the following equation:

{

2

Jˆ = min Φ − K M J J

+α J

2

}

(28)

where Φ is an N × 1 vector that contains the electrical potentials recorded at the scalp electrodes; J is a 3M × 1 vector representing the three orthogonal components of the current sources at M positions within the head conductor volume and KM is the field matrix for the realistic meshless model; α is the Tikhonov regularization parameter. The estimated source vector is then:

[

∧

]

J = TK M J = K M T K M K M T + αI + K M J = RJ

(29)

where R is the resolution matrix, which describes the deviation between real and estimated sources; ideally R should be the identity matrix. The basic idea of sLORETA consists on normalizing the estimation using block by block inverse of this solution matrix: ∧T

∧

J l (R ll )−1 J l

(30)

∧

where J l is a 3 × 1 vector containing the estimated source at the l-th voxel and containing the l-th diagonal of the solution matrix (Pascual-Marqui, 2002).

R ll

is a 3× 3 matrix

b) Signal Simulation The simulated EEG signals that were utilized for the tests were generated in a spherical head model through the algorithm presented by (Drozd et al., 2005); the signals are generated by a number of point sources (Si) placed inside the 3D homogeneous and isotropic conductive volume. The sources radiate their signals omnidirectionally. The signal sensors (electrodes) are placed in the volume surface, in a standard 10-20 configuration of 64 electrodes. The signals at the electrodes are lineal mixtures of the radiation from all the sources utilized. Each source generates a signal Si(n) defined by: S i (n) = − f i (n)ai (n), n = 0,K, N − 1

f i (n) =sin (2πn ⋅ 0.005 ⋅ (2000 − n ) / 2000 )

ai ( n ) =

1 2π

(

exp − (n − 500 )2 ⋅ 0.016

(31)

)

where fi(n) is a variable frequency signal amplitude modulated by a Gaussian bell ai(n) and N is the number of temporal samples of a trial. An isotropic, homogeneous, unitary radius spherical head conductor model is used. The summation of the signals coming from all sources to the sensors is carried out according to: L

Φ e ( n) =

∑ s ( n) d i

(32)

ei

i =1

181

where Φ e is the electrical potential in sensor (electrode), L is the number of signal sources (in our case L is equal to one). The distance coefficient dei for a pair consisting of the sensor at re and the source at ri is defined by:

⎧1, re − ri = 0 ⎪ d ei = ⎨ 1 ⎪ r − r , re − ri > 0 ⎩ e i

(33)

c) MRI Segmentation method. The anatomical information required for the construction of the realistic volume conductor model was obtained from MRI scans from the IBSR (Internet Brain Segmentation Repository), segmented through the mean shift procedure weighted by confidence maps, as reported by (Jiménez, et al., 2006). This segmentation procedure is a nonparametric estimation based on the mean shift (MS) algorithm that uses the local modes of the underlying joint space-range density function to define the voxel cluster centres. This scheme has proven to be robust under the presence of various types of noise at different levels and does not assume a particular probability density function for the voxel intensities. The quality of the class boundaries is improved by including an edge confidence map that represents the confidence of truly being in the presence of a border between adjacent regions. The confidence measure is also used to merge regions with weak edges, through the iterative application of transitive closure operations in a region adjacency graph (RAG). The automatic classification of the regions found by the MS filtering procedure, into white matter, CSF and gray matter regions is based on a digital brain atlas against which a spatial normalization is carried out. d) Localization Error Computation The proposed inverse solution was validated by estimating the localization differences ε which were computed by measuring the distance between estimated source location XM and simulated generators position XSPH:

ε = X SPH − X M

[

i i i XiSPH = xsph , ysph , z sph

[

XiM = xmi , ymi , zmi

]

]

(34)

2.3 Methodological Summary The methods and procedures used in this work (sections 2.1 and 2.2) can be summarized as follows: 1. 2. 3. 4. 5.

MRI segmentation to classify voxels and extract gray matter; this volume is used as a solution space for the inverse problem (section 2.2 c). EEG signal simulation generated by current sources placed at different positions; a 64 electrode set was defined, according to the 10-20 standard configuration (section 2.2 b). Computation of the lead field matrix KM for the realistic meshless model (Eq. 27) and of the field matrix K ∞ (Eq. 20) for an infinite homogeneous medium (section 2.1) Estimation of the inverse problem solution using sLORETA (section 2.2 a). Determination of localization errors (section 2.2d).

3. Results

182

Twenty five simulated EEG signals were generated following the procedure described in section 2.2b for several source locations inside the unitary spherical model. Brain anatomical information was obtained from IBSR, corresponding to a subject, without reported pathology, and registered to the Talairach space. To extract the gray matter volume, these MR images were segmented following the procedure described in section 2.2c; the voxels corresponding to gray matter ( 1 × 1 × 1 mm3 resolution) were subsampled in order to obtain an 7 × 7 × 8 mm3 voxel resolution; this produced 2472 gray matter voxels, with 854 of those belonging to the sphere. We assumed a homogeneous conductivity of 1.0 S/m2 and the volume was normalized, in order to obtain a maximum unitary dimension and to be able to register it with respect to the unitary surface, with which the simulated signals were generated. The center of the spherical volume conductor used for simulation was registered to the center of the Talairach space of the segmented gray matter, as can be seen in Fig. 1, where the gray matter voxels are represented by blue dots and the electrodes by red dots superimposed on the sphere.

Figure 1. Spherical volume conductor (used for the direct problem) registered with the center of the realistic gray matter conductor volume (used for inverse problem). Blue dots correspond to the gray matter surface; the 64 electrode location is also shown in red dots.

In order to find the most suitable parameters of the Gaussian weighting functions (R and λ) in Eq. 7, the localization errors were computed for different parameter combinations. This was carried out by simulating ten EEG signals, considering the source positions indicated in Table 1. The ROI radius (R) and the attenuation of the Gaussian kernel (λ) were varied according to Table 2, to obtain several meshless combinations (MC).

Table 1. Source positions used to simulate ten EEG signals (units relative to the unitary-radius sphere) Simulated signal

Position X

Position Y

Position Z

1

0.851

-0.080

0.253

2

-0.874

-0.172

0.253

3

-0.172

-0.161

0.805

4

-0.264

-0.770

0.069

5

0.103

0.793

0.069

6

0.483

0.471

0.345

7

0.287

-0.713

0.253

8

-0.632

0.563

0.253

9

0.425

0.655

0.253

10

-0.218

-0.448

0.437

183

Table 2. The different meshless combinations of parameters R and λ, used in tests. MC 1

MC 2

MC 3

MC 4

MC5

R

0.305

0.305

0.3529

0.3529

0.3529

λ

0.1411

0.1176

0.1411

0.1176

0.1626

R/ λ

2.17

2.6

2.5

3.0

2.17

Fig. 2 shows the normalized error computed for the ten simulated signals, after applying the tested parameter combinations for the weighting functions. In order to appreciate the variation of errors obtained for different source locations, we have plotted the linear trend for each parameter combination; it can be seen that for MC1 combination (R/λ =2.17) the error remains constant for all the tested positions. In all cases, sLORETA was applied using a regularization parameter α=0, to find the inverse solution.

70%

MC1 IHM

60% MC2 Localization Error

50%

MC3 MC4

40% MC5 30% 20% 10% 0% 1

2

3

4

5

6 EEG

7

8

9

10

Linear Trend IHM Linear Trend MC1 Linear Trend MC2 Linear Trend MC3 Linear Trend MC4 Linear Trend MC5

Figure 2. Localization error (ε) computed for the ten simulated signals, applying the four parameter combinations. The lines represent the linear trend of the error for each parameter combination and the infinite homogeneous medium (IHM).

In order to test the method’s performance against source depth, 12 signals were simulated gradually decreasing the distance between the source and the scalp or, equivalently, increasing the distance between the source and the space center (Normalized distance to center), as can be seen in Table 3.

184

Table 3. Source position used to simulate twelve EEG signals (units relative to the unitary-radius sphere). EEG

Normalized Distance to center

1

0.347

2

0.410

3

0.447

4

0.494

5

0.505

6

0.537

7

0.567

8

0.600

9

0.663

10

0.758

11

0.821

12

0.838

Fig. 3 shows the localization error obtained for these sources, plotted versus the distance from the sphere center or from the Talairach space center for the case of gray matter. The linear trend is also included, but only as an indication of the function error behavior vs. depth. This shows an apparent decrement of errors, as the source approaches the surface, whenever the parameter combination is R/λ=2.17 (MC1), while for the other combinations the error shows a more erratic behavior. Also, we observe that for the MC1 combination, the error stays below 30% for normalized distances above or equal to 0.4.

185

MC1

70%

IHM 60% MC2

Localization Error

50%

MC3 MC4

40% MC5 30%

Linear Trend MC1 Linear Trend IHM Linear Trend MC2 Linear Trend MC3 Linear Trend MC4 Linear Trend MC5

20%

10%

0% 0.3

0.4

0.5

0.6

0.7

0.8

0.9

Normalized Distance to Center

Figure 3. Localization error versus source depth (distance from space center to source).

Given the results shown in Fig. 2 and Fig. 3, in the last test a comparison was carried out only for the MC1 combination, against the behavior of source localization using an infinite homogeneous medium for fifteen additional EEG signals. Fig. 4 shows that localization errors of MC1 combination, in all cases, are not bigger than 30% (or 24 mm for our volume model); in contrast, the IHM presents, for two cases, localization errors bigger than 40%.

MC1

70%

Localization Error

60% 50%

IHM

40% 30%

Linear Trend MC1

20% 10% 0% 0.3

0.4

0.5

0.6

0.7

Normalized Distance to Center

0.8

Linear Trend IHM

Figure 4. Localization error vs source depth (distance from space center to source): comparison of MC1 combination and infinite homgeneous medium solution.

In order to test more extensively the performance of the proposed method, we simulated 200 different dipoles, located at randomly selected eccentricities. We obtained an average localization error as shown in Table 4. Table 4. Average localization error (mm) for 200 randomly located dipoles. Voxel resolution: 7×7×8mm3 Axis

Meshless Combination 1(R/ λ=2.17)

Infinite Homogeneous Medium

X

4.955 ± 3.477

3.91 ± 3.228

Y

10.775 ± 4.934

10.515 ± 4.805

Z

5.0 ± 4.351

4.72 ± 4.479

3D

13.959 ± 5.102

13.272 ± 5.041

186

Computation time for the proposed solution is dominated by the determination of the lead field matrix. In our simulations, this matrix is computed for 2472 nodes under Matlab, and it takes 16 minutes for the complete process in a standard Pentium IV machine running at 2.99 GHZ with 512 MB of system RAM.

4. Discussion A meshless method has been introduced to find a lead field matrix that characterizes an electromagnetic conductor medium (the head, in the case of EEG applications) to be used for inverse problem solutions. The proposed method has the main advantage of not needing the construction of a mesh that interconnects the nodes forming the volume; in order to calculate the potential on each node, it only requires that the quantity and position of the nodes in its neighborhood be determined (Eq. 8). The precision of the source localization is not affected and the complexity of the operations is considerably reduced. An advantage of the proposed methodology is that the solution space is always guaranteed to lie within the gray matter voxels, as the segmentation procedure used incorporates the priors from a digital brain atlas for anatomical structure classification. Throughout the evaluations that have been carried out, we have observed that the most relevant parameters that affect the method’s performance are the node’s region of influence, measured by its radius (R), and the attenuation (λ) of the contribution of the surrounding nodes. This led us to test different combinations of these two parameters (Table 2) for the construction of the lead field matrix to find the inverse solution; the tests were carried out by simulating several EEG signals (noiseless) under controlled, but random source positions (Table 1); the localization errors were computed following Eq. 34. Fig. 2 shows the trend of these errors for the ten simulated conditions, where it can be seen that the combination with the best performance is MC1 (R/λ=2.17), whose average error is 10%, in contrast with other combinations, where the error goes up to 40%. These results are consistent with those reported by (Ellenrieder et al., 2005) for the direct problem solution. It is clear that the meshless method needs a considerable quantity of nodes within a given ROI to achieve a good precision; so that with a greater region of influence (R), greater accuracy in the potential approximation is obtained, as can be observed in Fig. 2. In this sense, the tests with EEG signals for different depths of the current sources (Fig. 3), show that the performance of the meshless method with larger radius is better for the localization of deep sources; nonetheless, this precision is lost when the sources are located near the surface. This is verified when we compare the meshless combinations 1 and 3, where RMC3>RMC1. On the other hand, if we compare MC1 and MC2, which have same radii but different attenuations (λMC2 < λMC1), the localization error increases with a smaller attenuation, for source location near the surface. Also, parameter combinations having the same ratio (R/λ=2.17, Table 2) were compared, (MC1 and MC5), where RMC5>RMC1 and λMC1 < λMC5; the location error obtained with MC1 is lower, demonstrating that the parameter that has more impact in the results is the radius of the ROI (Figs. 2 and 3). Fig. 4 shows the performance for different source depths and allows the comparison of the selected MC1 combination versus the IHM model; it can be seen that the behavior of the two methods has a similar linear trend. MC1 presents less than 10 % of localization error, for sources near the surface, but for deeper sources the MC1 error gets closer to that of the IHM case. Computations for the infinite homogeneous model always show better localization performance, when compared to the meshless solution. This is actually expected from the fact that the simulated EEG signals were generated assuming an infinite homogeneous medium, and only the inversion procedure is introducing error, as compared to the local approximations performed in the meshless inversion. In Table 4, the average error obtained for 200 recordings at different depths can be observed. The average 3D error is 13 ± 5 mm, for a 7 × 7 × 8 mm3 resolution, which is similar to that reported in the literature for other approaches (Pascual-Marqui, 2002). It can also be appreciated that the error is not uniformly distributed, and that it has a higher magnitude along the y-axis. This is consistent with the fact that the virtual electrodes are located on a spherical geometry and thus the distance between the real anatomical model and the sphere would be larger around the temporal regions, which correspond to the y-axis direction. The volume conductor model proposed in this work is not the best possible approximation to a realistic brain model, which would consist on a multilayer, inhomogeneous one. Unfortunately, the potential approximations with shape functions that are being used are smooth, while in general the transitions of the normal gradient components at the layer interfaces are not. An alternative to be considered in future work for extensions to inhomogeneous media models might be on the line

187

suggested in (Ellenrieder et al., 2005), where in the forward formulation each layer is dealt with separately. Layer combination is achieved by establishing layer boundary conditions, at the expense of an increased size of the resulting linear system.

5. Conclusions A contribution of this work is the establishment of a sound methodology to solve the EEG inverse problem with a meshless method. Although this method uses a sparse lead field matrix, in contrast to the full matrix given by an infinite homogeneous medium solution, the performance is good for deep source localization, with the added benefit of integrating true anatomical information of the patient, which is not effectively carried out in the IHM method. The lead field matrix computation developed in this work, suggests an alternative construction of realistic volume conductor models for the EEG inverse problem solution that would eventually lead to better precision in the localization of electrical sources in ill-conditioned systems. An extension of this work would be to thoroughly analyze the regularization effects of the meshless model, since we have informally noted that, based on singular value analysis, the conditioning of KM was always better than that of K ∞ .

Acknowledgements We thank Dr. Joaquín Delgado from the Department of Mathematics at the Universidad Autónoma Metropolitana - Iztapalapa, for extensive reviews and helpful discussions regarding this paper.

References Akalm-Acar Z, Gencer N. An advanced boundary element method (BEM) implementation for the forward problem of electromagnetic source imaging. Physics in Medicine and Biology, 49: 5011-5028, 2004. Belytschko T, Krysl P, Krongauz Y. A three-dimensional explicit element-free galerkin method. International Journal for Numerical Methods of Fluids, 24:1253-1270, 1997. Cuffin N. A method for localizing EEG source in realistic head models. IEEE Transactions on Biomedical Engineering, 42(1): 68-71, 1995. Darvas F, Pantazis D. Mapping human brain function with MEG and EEG: methods and validation. Neuroimage, 23, S289-S299, 2004. de Munck J, Peters M. A fast method to compute the potential in the multisphere model. IEEE Transactions on Biomedical Engineering, 40(11): 1166-1174, 1993. Ding L, Lai Y, He B. Low resolution brain electromagnetic tomography in a realistic geometry head model: a simulation study. Physics in Medicine and Biology, 50: 45-56, 2005. Drozd M, Husar P, Nowakowski A, Henning G. Detecting evoked potentials with SVD- and ICA-based statistical models. IEEE Engineering in Medicine and Biology Magazine, 2005. Ellenrieder N, Muravchik C, Nehorai A. A meshless method for solving the EEG forward problem. IEEE Transactions on Biomedical Engineering, 52(2): 249-257, 2005. Ellenrieder N, Muravchik C, Nehorai A. Effects of geometric head model perturbations on the EEG forward and inverse problems. IEEE Transactions on Biomedical Engineering, 53(3): 421-429, 2006. Ellenrieder N, Muravchik C, Nehorai A. MEG forward problem formulation using equivalent surface current densities. IEEE Transactions on Biomedical Engineering, 52(7): 1210-1217, 2005. Eskola H. Utilization of MRI information in EEG studies. International Journal of Bioelectromagnetism, 1(1): 54-61, 1999. Hämäläinen M, Sarvas J. Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data. IEEE Transactions on Biomedical Engineering, 36(2): 165-171, 1989. He B, Musha T, Okamoto, Y, Homma S, Nakajima Y, Sato T. Electric dipole tracing in the brain by means of the boundary element method and its accuracy. IEEE Transactions on Biomedical Engineering, 34(6): 406- 414, 1987. Hédou-Rouillier V. A difference method to solve the forward problem in electroencephalography (EEG). Journal of Computational and Appied Mathematics, 167: 35-58, 2004. IBSR. Internet Brain Segementation Repository. http://www.cma.mgh.harvard.edu/ibsr. Jiménez J, Medina V, Yáñez O. Data-Driven brain MRI segmentation supported on edge confidence and a priori tissue information. IEEE Transactions on Medical Imaging, 25(1): 74-83, 2006. Li Z, Zhang Y, Zhu S, He B. Solving the ECG forward problem by means of a meshless finite element method. Physics in Medicine and Biology, 52: 287-296, 2007.

188

Mirkovic N, Adjouadi M, Yaylali I, Jayakar P. 3-D localization of epileptic foci integrating EEG and MRI data. Brain Topography, 16(2): 111-119, 2003. Mosher J, Leahy R, Lewis P. EEG and MEG: forward solutions for inverse methods. IEEE Transactions on Biomedical Engineering, 40(3): 245-259, 1999. Neilson L, Kovalyov M, Koles Z. A computationally efficient method for accurately solving the EEG forward problem in a finely discretized head model. Clinical Neurophysiology, 116: 2302-2314, 2005. Pascual-Marqui R, Esslen M. Functional imaging with low resolution brain electromagnetic tomography (LORETA): a review. Methods and Findings in Experimental and Clinical Pharmacology, Art. 24 C: 91-95, 2002. Pascual-Marqui R, Michel C. Low resolution electromagnetic tomography: A new method for localizing electrical activity in the brain. International Journal of Psychophysiology, 18: 49-65, 1994. Pascual-Marqui R. Review of methods for solving the EEG inverse problem. Internacional Journal of Bioelectromagnetism, 1(1): 75-86,1999. Pascual-Marqui R. Standardized low resolution brain electromagnetic tomography (sLORETA): technical details. Methods & Findings in Experimental & Clinical Pharmacology, 24D:5-12, 2002, Author’s version. Rao L, He R, Li Y, Liu S, Yan W. Computations of electroencephalography and magnetoencephalography for real head model. IEEE Transactions on Magnetics, 36(4): 1812-1816, 2000. Spinelli L, Gonzalez Andino S, Lantz G, Seeck M, Michel C. Electromagnetic inverse solutions in anatomically constrained spherical head models. Brain Topography, 13(2): 115-125, 2000. Srinivasan R. Methods to improve the spatial resolution of EEG. International Journal of Bioelectromagnestism, 1(1): 102-111, 1999. Teubner M, Nixon J, Rasser P, Bottema M, Clark R. Source localisation in a real human head. Brain Topography, 17(4): 197-205, 2005. Tissari S, Rahola J. Error analysis of a galerkin method to solve the forward problem in MEG using the boundary element method. Computer Methods and Programs in Biomedicine, 00: 1-14, 2003. Tripp J. “Physical concepts and mathematical models,” Biomagnetism: An interdisciplinary approach, S.J. Williamson, G.L. Romani, L. Kaufman, and I. Modena, Eds. New York: Plenum, 101-139, 1983. van ‘t Ent D, Munck J, Kaas A. A fast method to drive realistic BEM models for E/MEG source reconstruction. IEEE Transactions on Biomedical Engineering, 48(12): 1434-1443, 2001. Wu Q, Yan W, Shuen X, Li Y, He S, Xu G. EEG source reconstruction based on the boundary-element method and weigthted minimum norm approaches. IEEE Transactions on Magnetics, 39(3): 1547-1550, 2003. Xuan L, Zeng Z, Shanker B, Udpa L. Element-free galerkin method for static and quasi-static electromagnetic field computation. IEEE Transaction on Magnetics, 40(1): 12-20, 2004. Zhai Y, Yao D. A radial-basis function based surface laplacian estimate for a realistic head model. Brain Topography, 17(1): 5562, 2004. Zhang Y, Shao K, Xie D, Lavers J. Meshless method bases on orthogonal basis for computational electromagnetics. IEEE Transactions on Magnetics, 41(5): 1432-1435, 2005. Zhang X, van Drongelen W. High-resolution EEG: Cortical potencial imaging of interictal spikes. Clinical Neurophysiology, 114:1963-1973, 2003. Zou J, Xie Y, Yuan J, Cui M, Chen S. Source region contracting method for EEG source reconstructions. IEEE Transactions on Magnetics, 40(2): 1128-1131, 2004.

189