Solving the electroencephalography forward problem with a meshless method

N. von Ellenrieder¹ and C. H. Muravchik²

Laboratorio de Electrónica Industrial, Control e Instrumentación (LEICI),
Dto. Electrotecnia, Facultad de Ingeniería, U.N.L.P. CC 91, 1900 La Plata, Argentina.
¹ ellenrie@ing.unlp.edu.ar
² carlosm@ing.unlp.edu.ar

Abstract ¾ We describe a numerical method to solve the quasistatic Maxwell equations to obtain the electric potential distribution generated by a point source of current density inside a body of arbitrary shape and constant conductivity. The method needs only a set of nodes on the surface and inside the body, but it does not need a mesh connecting the nodes. The proposed meshless method is compared against the boundary elements method evaluating its performance when solving the electroencephalography forward problem.

Keywords ¾ Electroencephalography. Numerical Solution of Partial Differential Equations. Meshless Methods.

I. INTRODUCTION

One of the aims of electroencephalography (EEG) is to determine the location, orientation and intensity of the electrical activity generated by sources of neuronal activity in the brain. These sources can often be described as current density sources which generate an electric potential distribution on the surface of the head that can be measured with a suitable array of electrodes. The computation of the electric potential distribution based on sources of known location, orientation and intensity is known as the EEG forward problem. The estimation of source parameters based on electric potential measurements is the so-called EEG inverse problem. This paper essentially deals with the forward problem.

To get a good solution of the EEG forward problem it is necessary to correctly model the shape of the head layers and their respective conductivities (Hämäläinen and Sarvas, 1989). Due to its irregular geometry, this means that the problem has to be solved numerically. The most common numerical methods, such as finite elements (FEM) or boundary elements (BEM) require placing some nodes on the surface of the head layers or inside them, partitioning the domain of the problem into a set of small elements. A simple variation is assumed for the electric potential over each element, and the differential or integral equations are solved to get the potential at the nodes (de Munck, 1992; van Oosterom and Strackee, 1983).

While it is relatively easy to select the location of the nodes, based for example on magnetic resonance scans or in X-ray tomography, the generation of a mesh connecting them to define the elements is a difficult, time-consuming job, even if done automatically. This is especially true if the elements need some regularity, as is the case, for example, to keep a reasonable shape factor of the surface mesh for BEM (de Munck, 1992). To overcome this complications a number of methods that do not require a mesh connecting the nodes are being developed for a wide range of applications (Belytschko et al., 1997; Mukherjee and Mukerjee, 1997; Nayroles et al., 1992; Zhu et al., 1998).

In this paper we develop a meshless method for solving the forward problem of the quasi-static Maxwell equations. In particular, we will focus on a meshless method for the EEG forward problem utilizing a Local Boundary Integral Equation (LBIE) method (Zhu et al., 1998). Since at this point we are mainly interested in the viability of a meshless method for solving the EEG forward problem, we chose a simple model for the head, consisting in a constant conductivity arbitrarily shaped solid body. No detailed analysis of the errors and the numerical behavior of the method is presented here. An idea on how to demonstrate consistence and convergence for meshless methods can be found in (Duarte and Oden, 1996; Wendland, 1999).

In the following section the LBIE method is extended to solve the Poisson equation in a bounded three dimensional (3D) domain, with general boundary conditions and unknown boundary, except for a set of points. In the third section we address some complications arising when this method is applied to the EEG forward problem and restate it to avoid them. The fourth section focuses in the solution of the Laplace equation with Neumann boundary conditions (known flux across the boundary), based on which the EEG forward problem solution can be found. In the fifth section we present results obtained with this method and we compare them with results using BEM. The results are compared for a spherical head model. This simple geometry was chosen because its analytic solution to the forward problem is known (Zhang, 1995), and this facilitates the analysis of the numerical performance of the method and the comparison with other methods. Finally, in the sixth section we draw some conclusions with regards to the performance of the method and state the next steps in this line of work.

II. LBIE METHOD IN 3D DOMAINS

We have to obtain the electric potential 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 et al., 1993). Let W_G be the region of space occupied by the body and ¶W_G its boundary (see Fig. 1).

Figure 1: Layout of the electromagnetic problem

The LBIE method places nodes in points x_i inside or on the surface of the domain W_G and associates an equation to each node for the electric potential in it. A linear system of equations is posed whose solution is a vector with the potential at the nodes. The equations relate the potential F of node i to an arbitrary region surrounding the node. It is possible then to write an expression for the electric potential in a node by using Green's identity

where and are sufficiently smooth functions.

Let F denote the electric potential and Y the so-called fundamental solution, i.e. the solution to

Then, in Green's identity we get

(1)

The regions W_L could have arbitrary shape. In particular, when they are spheres centered in a node, the surface integrals in (1) take a very simple form. In addition, it is possible to make by taking , where r₀ is the radius of the sphere W_L. For internal nodes, i.e. those within W_G, it is always possible to find a value of r₀ small enough to make sure that the sphere W_L is completely contained in the domain W_G. Then (1) for internal nodes becomes

(2)

For nodes on the surface of the body, the local region W_L may be defined as the intersection between a sphere and the domain W_G (see Fig. 1). The boundary of the local region is formed by two surfaces: . The surface S_f is part of a spherical shell and S_r is part of the real boundary ¶W_G of the problem.

If the radius of the sphere defining the local region is small, the unknown surface S_r can be taken as a circular piece of a plane normal to n, the normal vector of the real surface. Using this approximation S_f is half a spherical shell and the vector x_i - x, x Î S_r is normal to n. From (1) we get for surface nodes

(3)

To solve the integral equations (2) and (3), the electric potential is locally approximated in the neighborhood of any point x Î W_G assuming a simple variation

(4)

where p(x) is the basis of the approximating functions. In this work we use a complete monomial basis of order 2, therefore

The vector a is a vector containing the coefficients of the approximation. Due to the local character of the approximation this vector varies in space (a = a(x)). It can be obtained by minimizing the following functional

(5)

The weight functions w_j(x) are responsible for the local character of the approximation. They give more importance to the nodes near to x, decrease as the nodes become farther apart from x and are null for the farthest ones. In this way, only the n nearest nodes to the point x participate in the determination of the coefficients a(x). There exist several options, but the weight functions used in this work are the so-called radial Gaussian functions defined as

(6)

where l_j is a positive constant and R_j is the size of the support of the weight function, also known as the influence region of node j.

Equation (5) posses a moving least squares problem for a, with the following well known solution

(7)

where

The coefficients vector is . The minimization of the functional (5) minimizes the difference between the coefficients and the electric potential at the nodes F(x_i), i = 1, ... n. To get a least squares problem with a unique solution, the number of nodes in the influence region of the point x has to be larger than the number m of elements from the basis of the approximation p(x) (m = 10 for the proposed monomial basis of order 2).

Let then, from (4) and (7)

(8)

The functions Y_j(x) are called shape functions of node j, and it can be shown that

Notice that since Y_j(x_i) ¹ d_ij the method does not interpolate the value of the electric potential between the nodes as BEM or FEM do. This explains why the coefficients usually have a different value than the electric potential F^h(x) at the nodes. Figure 2 shows a one- dimensional representation of this situation.

Figure 2: Electric potential approximation.

If the approximation F(x) = F^h(x) is considered valid, from (8) and (2) it is possible to write the electric potential of internal nodes as

(9)

and the potential of surface nodes, from (3), as

(10)

III. THE FORWARD PROBLEM IN EEG

The EEG forward problem has Neumann boundary conditions; i.e. the normal flux component through the body surface is null. The group of active neurons, source of the electric activity, can be described in many cases as a point source of current density or an electric potential dipole (de Munck et al., 1988). Hence the problem, in local or differential form, is the following

where q and p are the dipole moment (intensity and orientation) and position respectively, with p Î W_G, p Ï ¶W_G. Although it is possible to solve the forward problem using this expressions in (9) and (10), this approach has some drawbacks.

It can be seen from the equations that the dipole must be in one of the local regions, i.e. . This means that they should cover all the body () and this is done by making the local regions as large as the distance between neighboring nodes. It will be seen that the performance of the method for such a situation is not good.

Another problem arises because the real electric potential has spatially rapid variations in the neighborhood of the dipole, with values up to ±¥. The approximation to the real potential made with the proposed set of basis functions is necessarily very poor.

There is a means of avoiding these complications restating the forward problem as follows. The electric potential can be found as the superposition of two terms F(x) = F_F(x) + F_N(x), where each term is given by

(11)

The term F_F(x) corresponds to the electric potential generated by a dipole source in an infinite homogeneous medium. This term and its normal flux can be computed analytically hence, in order to solve the complete forward problem, the only term that needs to be computed numerically in (11) is F_N. This term corresponds to a Laplace equation with Neumann boundary conditions whose solution is smoother than the solution for the original problem.

IV. IMPLEMENTATION

In this section we describe the implementation of the proposed meshless method for solving the Laplace equation in the domain W_G with Neumann boundary conditions.

Writing an equation for each node based on expressions (9) and (10), with null Laplacian and known normal flux through the boundary, we get the following linear system of equations

(12)

where

Y_ij = Y_j(x_i)

For internal nodes:

and for surface nodes:

(13)

If we consider that in the expression (4) for the electric potential the coefficients a(x) are constant in a neighborhood U of the point x and that the sphere W_L is contained in the neighborhood U, then the shape functions are given by

and the rows of matrix G by

where . The integral now can be easily solved analytically, with a change of coordinates of the functions of the basis to local spherical coordinates.

For internal nodes the surface is a spherical shell and the integral gives

For surface nodes is half a spherical shell and if the coordinate system is chosen so that the normal vector of the surface equals to the unit vector along the axis x₃, the integral results in

We can write Y^T = p^T(x)M and working in local coordinates for node i, we get then

Y^T(x_i) = [1 0 0 0 0 0 0 0 0 0] M_i

Since the surface nodes are the only known points of the surface, the normal component of the flux through the boundary is also only known at the surface nodes. Again if we chose a small radius for the local regions, the surface integral in (13) can be solved by assuming a simple variation of the flux over the surface S_r. A logic choice would be to have a linear variation, for the flux is the gradient of the potential for which a quadratic variation is assumed in the local regions. Then the integral in (13) can be solved analytically

Thus, the coefficients from (12) must be chosen as those which are solution to the system

(14)

with the rows of matrix H given by

and the elements of vector f:

with k = 0 for internal nodes and k = 1 for surface nodes.

The local character of the methods is evident because only a certain number of neighboring nodes affect the potential of each point. This is the reason why the matrix H is a sparse matrix. Moreover, note that the potential problems with Neumann boundary conditions do not have a unique solution but rather an infinite number of them that differ only in a constant value. Hence, the matrix H has a null singular value, corresponding to a constant eigenvector. This means that the problem has to be solved with the generalized inverse or with a technique such as deflation.

Once the coefficients are known, the electric potential approximation for any point x Î W_G is given by the expression (8). For instance, if we are interested in the electric potential at the nodes we have

where is a vector containing the computed potential at the nodes.

V. RESULTS

Several tests were setup to determine the performance of the proposed method. Although the method can be used for an arbitrarily shaped head model, in this section a spherical geometry was chosen, with radius and conductivity equal to one. As mentioned in the introduction, the reason for this geometry is that it allows for an analytical computation of the solution (Zhang, 1995).

For all the examples shown in this section, a Gaussian weight function (6) is used, with parameter set to R_j/l_j = 2.5. Another value for the parameter or different weight functions would result in different shape functions. Although a good choice of the weight function can improve the performance, this is not a crucial point. The shape function would still be smooth radial basis functions, with similar properties regarding consistence and convergence (Duarte and Oden, 1996).

The measure used for the electric potential distribution error on the surface of the sphere is the relative difference measure (RDM) (Schlitt et al., 1995)

where the sum runs over all the surface nodes, F_a is the analytically computed potential and F_n the numerically computed potential.

Figure 3: Convergence study - Tangential dipole

Figure 4: Convergence study - Radial dipole

The nodes where chosen in a roughly regular pattern over the surface and in the interior of the sphere, trying to maintain a constant distance between neighboring nodes. The rationale behind this distribution of nodes is to keep the uniformity of the predictable numeric errors (Wendland, 1999).

The convergence of the method was studied in first place; the results are shown in Figs. 3 and 4. The figures show the mean error versus the number of surface nodes of the problem for dipoles located at depths of 0.6 and 0.8 times the radius of the sphere. Figure 3 corresponds to the results obtained for tangentially oriented dipoles, and Fig. 4 to radially oriented ones. We performed a comparison against BEM, a method widely used to solve the EEG forward problem. The convergence of BEM for problems with the same number of surface nodes is also plotted in Figs. 3 and 4. It can be seen that for dipoles near the surface the performance of BEM is better, while for deeper dipoles both methods have a similar performance.

Figures 5 and 6 show the error as a function of the dipole depth. The problem was solved with 522 nodes, 162 of them on the surface; and for 1591 nodes, 642 of them on the surface. The error of the solution to the same problem solved with BEM is plotted too. It was computed for a sphere with 642 surface nodes and a mesh connecting them forming 1280 triangular elements over which the potential varies linearly. It is clear that BEM outperforms the proposed meshless method for dipoles near the surface, nevertheless the performance of the LBIE method is acceptable, and by increasing the number of nodes the error for dipoles near the surface can be reduced.

Figure 5: Error as a function of dipole depth. Tangential dipole.

Figure 6: Error as a function of dipole depth.
Radial dipole.

Figure 7: Sensitivity to the size of local regions.
r = radius of local region.

Figures 7 and 8 show the sensitivity of the method to the radius of the local regions and to the size of the influence region of the nodes. The latter is directly related to the number of nodes in the zone of influence. The results shown are from a sphere with 1591 nodes, 642 of them on the surface. In Fig. 7 we see that the sensitivity to the size of the local regions is very low, while Fig. 8 shows a large sensitivity to the size of the influence region. The size of the influence region was chosen as to include 50, 80 or 110 nodes, with fixed node density. If the influence region is small the method is too local, resulting in a poor global solution. On the other hand, if there the influence region is too big, the error of the approximating functions rises because the lineal approximation for the potential variation is inadequate.

Figure 8: Sensitivity to the size of influence regions.
N = number of nodes in influence regions.

CONCLUSIONS

In this work we described a numerical method for solving the EEG forward problems in bodies of constant conductivity of arbitrary shape. The method does not need a mesh connecting the nodes of the problem.

From the comparison against the boundary elements method, it is clear that the meshless method needs more nodes to achieve solutions with the same precision. This could have been predicted since BEM needs only surface nodes, roughly lowering the dimension of the problem in one order. On the other hand, the matrix related to the meshless method is sparse, while BEM matrices are full. This amounts to an important saving in time and in the number of operations needed to solve a system like (14). Since the inverse problem usually involves many consecutive computations of the direct problem, the savings in the computational effort becomes an important advantage. The FEM and meshless method matrices share the sparsity property therefore, a comparison between them would be interesting.

It is important to remember that the generation of the meshes for BEM and FEM in 3D problems is a computationally expensive job that is not required by the meshless method. This is especially true for the intricate shapes involved in the layers for the EEG and related problems.

Since the performance of the method is quite acceptable, the next step would be to extend it to problems with several layers of different conductivity. This is required if a better model of the head is used. An extension of the method to allow the computation of the magnetic field generated by the source is also planned, as a tool for solving the Magnetoencephalography (MEG) inverse problem.

Acknowledgement
This work was supported by Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Comisión de Investigaciones Científicas de la Pcia. de Buenos Aires (CICPBA) and the Universidad Nacional de La Plata; with grants from CONICET and Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT).

REFERENCES
1. Belytschko, T., P. Krysl and Y. Krongauz, "A Three-Dimensional Explicit Element-Free Galerkin Method", Int. J. Numer. Meth. Fluids 24, 1253-1270 (1997).         [ Links ]
2. de Munck, J. C., B. W. van Dijk and H. Spekreijse, "Mathematical Dipoles are adequate to describe Realistic Generators of Human Brain Activity", IEEE Trans. Biomed. Eng. 35-9, 960-966 (1988).         [ Links ]
3. de Munck, J. C., "A linear Discretization of the Volume Conductor Boundary Integral Equation Using Analytically Integrated Elements", IEEE Trans. Biomed. Eng. 39-9, 986-990 (1992).         [ Links ]
4. Duarte, C. A. M. and J. T. Oden, "Hp Clouds-an hp Meshless Method", Numerical Methods for Partial Differential Equations 12, 673-705 (1996).         [ Links ]
5. Hämäläinen, M. S. and J. Sarvas, "Realistic Conductivity Geometry Model of the Human Head for Interpretation of Neuromagnetic Data", IEEE Trans. Biomed. Eng. 36-2, 165-170 (1989).         [ Links ]
6. Hämäläinen, M. S., R. Hari, R. J. Ilmoniemi, J. Knuutila and O. V. Lounasmaa, "Magnetoencephalography - Theory, Instrumentation and Applications to Noninvasive Studies of the Working Human Brian", Reviews of Modern Physics, 65-2, 413-497, (1993).         [ Links ]
7. Mukherjee, Y. X. and S. Mukherjee, "The Boundary Node Method for Potential Problems", Int. J. Numer. Meth. Eng. 40, 797-815 (1997).         [ Links ]
8. Nayroles, B., G. Touzot and P. Villon, "Generalizing the finite element method: Diffuse aproximation and diffuse elements", Computational Mechanics 10, 307-318 (1992).         [ Links ]
9. Schlitt, H. A., L. Heller, R. Aaron, E. Best, and D. M. Ranken, "Evaluation of Boundary Element Methods for the EEG Forward Problem: Effect of Linear Interpolation", IEEE Trans. Biomed. Eng. 42-1, 52-58 (1995).         [ Links ]
10. van Oosterom, A. and J. Strackee, "The Solid Angle of a Plane Triangle", IEEE Trans. Biomed. Eng. 30-2, 125-126 (1983).         [ Links ]
11. Wendland, H., "Meshless Galerkin methods using radial basis functions", Mathematics of Computation 68, 1521-1531 (1999).         [ Links ]
12. Zhang, Z., "A Fast Method to Compute Surface Potentials Generated by Dipoles within Multilayer Anisotropic Spheres", Physics in Medicine and Biololy 40, 335-349 (1995).         [ Links ]
13. Zhu, T. , J. D. Zhang and S. N. Atluri, "A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach", Computational Mechanics 21 , 223-235 (1998).         [ Links ]