The EEG forward problem: theoretical and numerical aspects

D. Rubio¹ and M. I. Troparevsky²

¹ Escuela de Ciencia y Tecnología, Univ. Nac. de General San Martín, Argentina.
Instituto de Ciencias, Univ. Nac. de Gral. Sarmiento, Argentina. arubio@ungs.edu.ar
² Facultad de Ingeniería, Univ. de Buenos Aires, Argentina. mitropa@fi.uba.ar

Abstract — We focus on the Forward Problem of electroencephalography, discuss a mathematical model and state properties of its weak solutions. A static and a time-dependent model for the source are considered. Numerical solutions, obtained by a Boundary Element Method technique, are compared with the analytical ones and with EEG recordings.

Keywords — EEG. Forward Problem. Epilepsy. Modeling.

I. INTRODUCTION

The electrical activity inside the brain consists of currents generated by biochemical sources at cellular level. This activity can be measured by an electroencephalograph. In the case of epilepsy there are small zones inside the brain that give major contribution in the generation of the electric field. Neurologists have been interested in determining the location of the epileptogenic zones from the measured potential on the scalp in order to avoid invasive techniques.

The problem of determining the source location from the EEG waves is known as the Inverse Problem of Electroencephalography (EEG). In order to solve the inverse problem, we need to have an appropiate model of the Forward Problem (FP) of EEG that consists of calculating the EEG signal when we know the source location.

A typical mathematical model that describes this process is a PDE-boundary value problem of second order, based on the static approximation of Maxwell's Equations (see Hamalainen et al., 1993).

In this work we consider a simplified model of the human head (spherical). We calculate approximated solutions and compare them with the analytical one that exists in the spherical case.

Regarding the source, dipole source models as well as spatially distributed models can be found in Schimpf et al. (2002), Lagerlund (1999), de Munck (2002), Yetik et al. (2005a) and Yetik et al. (2005b). One dimensional source distribution to model the primary cortical response to nerve stimulation is theoretically analyzed in Nolte and Curio (2000). A spatiotemporal source analysis based on the spatiotemporal noise covariance matrix is developed in Huizenga et al. (2002). We also propose a time-dependent model for the electric source that approximates a dipole and compare simulated results with real data obtained from EEG recordings provided with by Centro Municipal de Epilepsia, Hospital Ramos Meja, Buenos Aires, ARGENTINA.

The paper is organized as follows. In section II the PDE mathematical model is presented. Weak solutions of the FP and their properties are introduced in section III. Section IV includes derivation of the integral equations and its discretization. Different models for the electric source are also described. In section V we present numerical simulations. The approximated solutions for the static dipole model are compared with the analytical ones. The simulated results in the case of a time-dependent source modelare compared with real data from EEG recordings. EEG signals plots and 3D plots illustrate the results.

II. THE MATHEMATICAL MODEL

The electric sources inside the brain produce electric and magnetic fields that can be modeled by the Maxwell Equations (see Hamalainen et al., 1993).

The electric current is assumed to be of the form

J(x)= σ(x) E(x) + J_i(x), (1)

where σ is the conductivity function, σE is the macroscopic electric field and J_i is the impressed current (microscopic level). During an epileptic seizure, spikes can be observed along the EEG signals. They are mainly produced by the impressed current.

Due to the high speed of propagation of the electric waves inside the head, there is no delay in the data captured by the EEG recorder. Hence, in order to find the location of the impressed current, we consider the EEG data at the instant at which one of the spikes achieves its highest amplitude. Therefore, a time-independent Maxwell equation may be used to model the relationship between the electric potential u and the impressed current J_i (see Hamalainen et al., 1993), that is,

∇ · (σ(x)∇u(x)) = ∇ · J_i(x) x ∈ G (2)

(3)

subject to

(4)

where G represents the head having different compartments with transition surfaces named S_i and [·] denotes the difference between the values of the functions inside the brackets through the indicated surface.

A. Models for the Head and the Conductivity Function

The model of the human head consists in three homogeneous, nested sets denoted by G₁,G₂,G₃, representing the brain, the skull and the scalp, G = G_j. The surface of transition between G_i and G_i+1, i =1, 2, is denoted by S_i. The conductivity function is approximated by a piecewise constant function of the form

(5)

where σ₁, σ₂, σ₃ are constant positive values.

For this model there exists an analytic solution when G_i are nested spherical volumes (Zhang, 1995; Mosher et al., 1999).

III. WEAK SOLUTIONS OF THE FORWARD PROBLEM

In this section we state the weak formulation of Eq. (2) and some properties of its weak solutions.

Suppose that u is a solution of Eq. (2) and multiply this equation by a function v. Assuming that both, u and v are regular enough to apply integral theorems to the resulting equation, u must verify the following identity

(6)

or, equivalently,

(7)

Note that the above integral equation requires only derivatives of first order while the differential equation (2) requires derivatives of second order.

A weak solution of Eq. (2) with boundary condition (3) is a function u that verifies Eq. (6) for all functions v with weak derivative of first order, i.e., v ∈ H¹(G), where

H¹(G)= { v ∈ L²(G) /∃w ∈ L²(G) with

(8)

with is the set infinitely differentiable functions with compact support in G. The function w that satifies (8) is called the generalized derivative of first order of v.

Remark III..1 Note that for the model described above, Eq. (7) is equivalent to

Since any solution of Eq. (2) is a weak solution we study some properties of the weak solutions referred to the model we have adopted.

Note that since the head is not isotropic, nor homogeneous media, an error is introduced when a piecewise constant function like the one defined in Eq. (5) is considered to model the conductivity function.

In Troparevsky and Rubio (2003a) we present the existence and uniqueness of the weak solutions and we proved the continuity of the weak solution with respect to the domain.

Theoretical and numerical results about the continuity of the weak solution with respect to the conductivity functions can be found in Troparevsky and Rubio (2003b).

IV. NUMERICAL APPROACH

A 3D Finite Element Method (FEM) can be chosen to solve Eq. (2) in the volume G, by using the weak formulation given in Eq. (6) (see for example, Bradley et al., 2001). Another approach frequently used to solve this problem is a Boundary Element Method (BEM)(see Brebbia et al., 1984). It requires to transform the equation defined in G into surface integrals, reducing the problem from 3D to 2D (see for example, Meijs et al., 1989; Schlitt et al., 1995; Ferguson and Stroink, 1997). Meshless methods were also developed for solving the Forward Problem of EEG (see von Ellenrieder et al., 2005).

Here we choose a BE Method that generates a full but much smaller matrix than the one that FEM techniques do.

A. From Volumes to Surfaces

Let u be the solution of Eq. (2) with boundary and interface conditions (3) and (4) where G and the conductivity function are defined as in the preceding section.

Let v be the solution of

Δv = δ(r₀) with r₀ ∈ G,

where Δ is the Laplacian operator and δ is the delta-Dirac distribution.

Multiplying Eq. (2) by v(r) we can derive the following Integral Equation for the potential u (see Hamalainen et al., 1993 and Sarvas, 1987),

(9)

for , where

(10)

Note that the only term containing information of the impressed current J_i is u₀(r), on the right hand side.

As we need to know the value of the potential u only on the external surface, we let r → r_j ∈ S_j obtaining an integral equation over the transition surfaces S_j ,

(11)

B. Discretization

There are many different techniques to discretize the integral equation (11). In Schiltt et al. (1995) and Meijs et al. (1989), for instance, a polyhedral surface Ŝ_j, composed by plane triangles is chosen to approximate the surfaces S_j. With respect to the discretization of the potential u over the chosen grid, in Schiltt et al. (1995) and Bradley et al. (2001), there are several comparisons of the results obtained for different approximations (see the references cited therein).

We discretize the integrals dividing each surface S_j into small spherical elements E_kj, i.e., S_j = E_kj . In addition, on each element E_kj the electric potential is approximated by the average of the values at the vertices of the element (see Rubio and Troparevsky, 2003) yielding

(12)

where C(u, j, k) is a constant that depends on the values of u at the nodes of E_kj.

This type of equation is obtained from different approximations like the COG method, the vertex method on triangular elements (see Schiltt et al., 1995) and by second-order interpolation functions (see Frinjs, 2000).

The integrals on the left hand side of (12) are the solid angle subtended by E_kj at r (see Santaló, 1973). These integrals become improper when r ∈ E_kj. To approximate its values, we use a geometric property that states that for any closed and smooth surface S, the solid angle at any of its points is -2π. We first approximate the integral on E_kj for r ∉ E_kj . Then we apportion the remaining value

among the elements E_kj that satisfy r ∈ E_kj , proportionally to the their areas (see Rubio and Troparevsky, 2004).

The last term to be discretized is the right hand side of Eq. (11) defined in (10), which requires a model for the impressed current J_i. We will discuss it in the next subsection.

The final discretized system takes the form

(D - A) · u = C (13)

where D is the diagonal matrix arising from the right hand side of (12), of the form

(14)

where , for i = 1, 2, 3, I is the identity matrix, A is the matrix resulting from the discretization of the surface integrals appearing in (12) and C is the discretization of u₀(r) (10).

Note that the values of u at the nodes of the three surfaces are involved in the calculation of u at any nodal point.

C. Impressed Current

Some neurologists state that in the case of epilepsy, the primary current J_i is concentrated in a small area around the source location, namely r_q.

Different models for the electric source may be found in the literature. Dipole source models and localization methods can be found in Schimpf et al.(2002) and de Munck et al. (2002). In Yetik et al. (2005a) and Yetik et al. (2005b) possible models for a stationary source are discussed and statistical selection methods are used to characterize them. One dimensional source distribution to model the primary cortical response to nerve stimulation is theoretically analyzed in Nolte and Curio (2001). In Huizenga et al. (2002) a spatiotemporal source analysis based on the spatiotemporal noise covariance matrix is developed. Estimation of stationary dipoles from MEG and EEG noisy data is discussed in de Munck et al. (2002). Other references can be found in the cited papers.

We propose two models for the impressed current. The first one is a static model that approximates a dipole, i.e., it is nearly zero outside a small area around rq

J_i(r) = q(r)δ(r - r_q) r ∈ G. (15)

The other model intends to reflect the electrical activity of the brain for a short period of time T around a spike of a seizure. We consider a function that depends on time, t, and position, r, and approximates a dipole (see Rubio and Troparevsky, 2005),

(16)

where

with α₁ = 100, α₂ = 3000, α₃ = 10000, α₄ =0.5, β₁ = 0.001, β₂ = 2500, β₃ = 2000, β₄ = 100 and q a fixed orientation.

With these source models we solve the FP for some instants within T trying to reproduce the seizure.

V. NUMERICAL SIMULATIONS

We simulate epileptic seizures by solving the discrete system (13) on the spherical domain. We place the origin of coordinates in the center of the sphere. The radii of G_i, i =1, 2, 3 are taken to be 0.071m, 0.078m and 0.085m. The conductity values in are σ₁ = 0.33, σ₂ = 0.0042 and σ₄ =0.33 (see Geddes and Baker, 1967).

For the case of the static model, the solution of Eq. (13) using a primary current of the form (15) is compared with the analytic solution for different source locations r_q and orientation q.

The solution obtained for the time-dependent model (16) is compared with real data obtained from the recordings of spontaneous activity during a seizure (provided by Centro Municipal de Epilepsia, Hospital Ramos Mejía, Buenos Aires, ARGENTINA).

A. Simulated vs. Analytic Solution

The solution for spherical domains can be calculated analytically (see Zhang, 1995 and Mosher et al., 1999), we compare its values on the scalp with the ones obtained by BEM.

We numerically solve Eq. (12) with the impressed current modeled by (15). Two quantitive measures for the numerical errors of the simulated scalp potential, are calculated; the relative difference measure (RDM) and the normalized relative difference measure (RDM*, also called NRDM), defined by

	(17)
	(18)

(see Meijs et al., 1989, von Ellenrieder et al., 2005). Taking 760 nodes per surface, for the source location = (0., 0.0062, 0.04) with orientation q₁ = (1.2, 0.6, 0.6)* 10^-9, we have RDM = 0.214602, RDM* = 0.206271, and for = (1, 0, 0)* 10^-5 with orientation q₂ = (1, 0, 0) * 10^-8, we obtain RDM = 0.0254408, RDM* = 0.00580414. We want to point out that these error measures decrease with the number of nodes.

B. Simulated Potential vs. Real Data

EEG signals taken from patients and their medical diagnosis were provided by neurologists from Centro Municipal de Epilepsia, Hospital Ramos Mejía, Buenos Aires. We consider the EEG waves sampled at 200 Hz. shown in Fig. 1. They correspond to a short period of time around a spike. Based on the information available we decided the source location r_q to be r_q = (1.2, 0.6, 0.6)cm and the source orientation q = 10^-10(1.2, 0.6, 0.6)nAm.

Figure 1: EEG Recordings.

A time period of 0.2 seconds centered in a spike (boxed in Fig. 1) was considered. We simulate the scalp electric potential with BEM by solving the discretized system (13) where the electric source was considered to follow the time-dependent law (16). This means that, at each instant t, we evaluate (16) and the result is considered as the source in (10) for the Eq. (13). Then we plot the potential distribution on the scalp for both, the real and the simulated case. With the latter we generate EEG signals and we also plot the resulting values on a sphere for fixed time instants.

EEG Signal Plots

In this subsection we present the EEG signals (real one and the simulated ones), i.e., the values of the potential u on each electrode channel varying in time (see Fig. 2). We simulate the electric potential, consider the electrode positions and plot the values of the potential at each channel (see Fig. 3).

Figure 2: The real EEG signal.

Figure 3: The simulated EEG signal.

We point out that the values of the numerical potential at the electrode positions are obtained by interpolating the simulated values, since we numerically solve the equation for potential at the grid nodes, which do not coincide with the electrode positions. Although there are some differences between the real and the simulated data, the shapes of the signals look similar around the spike for most of the channels.

3D Plots of the Potential

Another way of comparing the real vs. the simulated data is to plot the scalp potential distribution on a sphere at any fixed time instant. In Fig. 4 we present the 3D plot of the real EEG data at the spike time marked in Fig. 1. The approximated potential values

Figure 4: Measured Potential at a Spike Instant.

at that instant is plotted in Fig. 5. The color scale refers to the intensity of the electric scalp potential. We recall that real and numerical potential values

Figure 5: Simulated Potential at a Spike Instant.

are given at electrode positions and grid nodes, respectively. We use MATLAB interpolation to plot them on the spherical surface.

Some differences can be observed between the real and the simulated potential plots, however the pattern is similar. Moreover, the sense and direction of the gradient of the potential coincide. We will use the source model (16) with the location and strength of the electric source as an initial guess for the Inverse Problem in future works related to this patient.

VI. CONCLUSIONS

This work presents a time-dependent source model that seems to be suitable for the real case presented. This model can be adapted to different seizure processes by adjusting the parameters q, r_q, the constants α_j, β_j and by allowing r_q to change in time.

In order to calculate the approximated potential distribution we use estimations of the position r_q and the orientation q of the source based on clinical information. We approximate the EEG signal and plot it on a spherical head model. Then we compare the real and simulated data as EEG waves and as 3D plots. The results obtained with this source model encourages us to use it in future works to estimate the parameters of the model by means of an Inverse Problem Technique.

We point out that although the source model (16) depends on time, we use a FP static model, hence the dynamics of the problem is only partially shown.

VII. ACKNOWLEDGEMENTS
This work was partially supported by Fundación Antorchas under Grant 13900-2. The authors also want to thank the Centro Municipal de Epilepsia of Hospital Ramos Mejía, Buenos Aires, Argentina, for their collaboration and encouragement.

REFERENCES
1. Bradley, C.P., G. M. Harris and A. J. Pullan, "The Computational Performance of a High-Order Coupled FEM/BEM Procedure in Electropotential Problems", IEEE Trans. Biomedical Engineering, 48, 1238-1250 (2001).         [ Links ]
2. Brebbia, C.A., J.C.F. Telles and L.C. Wrobel, Boundary Element Techniques: Theory and Applications in Engineeering, Springer Verlag, N.Y. (1984).         [ Links ]
3. de Munck, J.C., H. Huizenga, L.J. Waldrop and R.M. Heethaar, "Estimating Stationary Dipoles from MEG/EEG Data Contaminated with Spatially and Temporal Correlated Background Noise", IEEE Trans. On Signal Processing, 50, 1565-1572 (2002).         [ Links ]
4. Ferguson, A.S. and G. Stroink, "Factors affecting the accuracy of the Boundary Element Method in the Forward Problem I: Calculating Surface Potentials", IEEE Transactions on Biomedical Engeneeering, 44, 1139-1155 (1997).         [ Links ]
5. Frinjs, J.H.M., "Improving the Accuracy of the Boundary Element Method by the Use of Second-Order Interpolation Functions", IEEE Transactions on Biomedical Engineering, 47, 1336-1346 (2000).         [ Links ]
6. Geddes, L.A. and L.E. Baker, "The specific resistant of biological materials-a compedium of data for the biomedical engineering and physiologist", Med. Biol. Eng. (1967).         [ Links ]
7. Hamalainen, M., R. Hari, R.J. Ilmoniemi, J. Knuutila and O. Lounasmaa, "Magnetoencephalography, Theory, Instrumentation and Applications to Noninvasive Studies of the Working Human Brain", Reviews of modern Physics, 65, 414-487 (1993).         [ Links ]
8. Huizenga, H., J.C. de Munck, L.J. Waldrop and R.P.P.P. Grasman, "Spatiotemporal EEG/MEG Source Analysis Based on a Parametric Noise Covariance Model", IEEE Trans. Biomedical Engineering, 49, 533-539 (2002).         [ Links ]
9. Lagerlund, T.D., "EEG Source Localization (Model-Dependent and Model-Independent Methods), Electroencephalography: Basic Principles, Clinical Applications, and Related Fields", Chap. 46, 1999, pp 809-822.         [ Links ]
10. Meijs, J.W.H., O.W. Weier, M.J. Peters, and A. Van Oosterom, "On the Numerical Accuracy of the Boundary Element Method", IEEE Trans. Biomedical Engineering, 36, 1038-1049 (1989).         [ Links ]
11. Mosher, J.C., R.M. Leahy and P.S. Lewis, "EEG and MEG: Forward Solutions for Inverse Methods", Trans. Biomedical Engineering, 46, 245-259 (1999).         [ Links ]
12. Nolte, G. and G. Curio, "Current Multuipole Expansion to Estimate Lateral Extent of Neuronal Activity; A theoretical Analysis", IEEE Trans. Biomedical Engineering, 47, 1347-1355 (2000).         [ Links ]
13. Rubio, D. and M.I. Troparevsky, "On the Accuracy of an Approximated Solution of the EEG Forward Problem", Proceeding of the RPIC 2003, San Nicolás, ARGENTINA, 1, 19-23, (2003).         [ Links ]
14. Rubio, D. and M.I. Troparevsky, "On the Approximation of the Auto Solid Angle for Solving Integral Equations", WSEAS Transactions on Mathematics, Issue 1, 3 (2004).         [ Links ]
15. Rubio, D. and M.I. Troparevsky, "Time Dependent Source Models for EEG", WSEAS Transactions on Mathematics (to be published) (2005).         [ Links ]
16. Santaló, L.A., Vectores y Tensores con sus Aplicaciones, EUDEBA, B. A., Argentina (1973).         [ Links ]
17. Sarvas, J., "Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem", Phy. Med Biol., 32, 11-22 (1987).         [ Links ]
18. Schlitt, H., L. Heller, R. Aaron, E. Best, M. Ranken, "Evaluation of Boundary Element Methods for the EEG Forward Problem: Effect of Linear Interpolation", IEEE Trans. Biomedical Engineering, 42, 52-58 (1995).         [ Links ]
19. Schimpf, P.H., C. Ramon and J. Haneisen, "Dipole Models for EEG and MEG", IEEE Trans. Biomedical Engineering, 49, 409-418 (2002).         [ Links ]
20. Troparevsky, M.I. and D. Rubio,"On the Weak Solutions of the Forward Problem in EEG", Journal of Applied Mathematics, 12, 647-656 (2003a).         [ Links ]
21. Troparevsky, M.I. and D. Rubio, "Weak Solutions of the Forward Problem in EEG for Different Conductivity Values", Journal of Mathematical and Computer Modelling (to be published) (2003b).         [ Links ]
22. von Ellenrieder, N., C.H. Muravchik and A. Nehorai, "A Meshless Method for Solving the EEG Forward Problem", IEEE Trans. Biomedical Engineering, 52, 249-257 (2005).         [ Links ]
23. Yetik, I.S., A. Nehorai, J. Lewinw and C.H. Muravchik, "Distinguishing Between Moving and Statiopnary Sources Using EEG/MEG Measurements with Application to Epilepsy", IEEE Trans. Biomedical Engineering (to be published) (2005a).         [ Links ]
24. Yetik, I.S., A. Nehorai, C.H. Muravchik and J. Hauuesen, "Line-Source Modeling and Estimation with Magnetoencephalography", IEEE Trans. Biomedical Engineering (to be published) (2005b).         [ Links ]
25. Zhang, Z., "A Fast Method to Compute Surface Potentials Generated by Dipoles within Multilayer Anisotropic Spheres", Phys. Med. Biol., 40, 335-349 (1995).         [ Links ]

Received: September 21, 2005.
Accepted for publication: February 2, 2006.
Recommended by Guest Editors C. De Angelo, J. Figueroa, G. García and J. Solsona.