## Servicios Personalizados

## Articulo

## Indicadores

- Citado por SciELO

## Links relacionados

- Similares en SciELO

## Compartir

## Latin American applied research

##
*versión impresa* ISSN 0327-0793

### Lat. Am. appl. res. v.36 n.2 Bahía Blanca abr./jun. 2006

**The EEG forward problem: theoretical and numerical aspects**

**D. Rubio ^{1} and M. I. Troparevsky^{2}**

^{1} *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*

^{2}

*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*_{1},*G*_{2},*G*_{3}, representing the brain, the skull and the scalp, G = *G _{j}*. The surface of transition between

*G*and

_{i}*G*

_{i}_{+1},

*i*=1, 2, is denoted by

*S*. The conductivity function is approximated by a piecewise constant function of the form

_{i}(5) |

where σ_{1}, σ_{2}, σ_{3} 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*^{1}(G), where

*H*^{1}(*G*)= { *v* ∈ *L*^{2}(*G*) /∃w ∈ *L*^{2}(*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*_{0}) with *r*_{0} ∈ *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*

_{0}(

*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*obtaining an integral equation over the transition surfaces

_{j}*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*. With respect to the discretization of the potential u over the chosen grid, in Schiltt

_{j}*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*. In addition, on each element

_{kj}*E*the electric potential is approximated by the average of the values at the vertices of the element (see Rubio and Troparevsky, 2003) yielding

_{kj}(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*. 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

_{kj}*E*for

_{kj}*r*∉

*E*. Then we apportion the remaining value

_{kj}

among the elements *E _{kj}* that satisfy

*r*∈

*E*, proportionally to the their areas (see Rubio and Troparevsky, 2004).

_{kj} 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*_{0}(*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 α_{1} = 100, α_{2} = 3000, α_{3} = 10000, α_{4} =0.5, β_{1} = 0.001, β_{2} = 2500, β_{3} = 2000, β_{4} = 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 σ

_{1}= 0.33, σ

_{2}= 0.0042 and σ

_{4}=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} = (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*_{2} = (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*= (1.2, 0.6, 0.6)

_{q}*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*to change in time.

_{q} 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.**