Revista Mexicana de Ingeniería Biomédica

Introduction

Different non-invasive techniques for brain scans have been developed using mathematical models. These include positron emission tomography, magnetic resonance imaging and electroencephalography. This final technique is the present study focus. Electroencephalography is the best known among the non-invasive brain investigation methods. It is based on the use of scalp located electrodes to record brain electrical activity. This recording is known as the Electroencephalogram (EEG). The electrical activity is generated by the bioelectrical activities of large neuron populations working synchronously [1, 2]. The EEG technique allows us to detect anomalies in the brain (damage, malfunction, etc.) which have traditionally been done by different diagnostic techniques. The problem is studied through a boundary value problem, which is obtained using a model that describes the head as a conductive layer medium. This model allows finding relationships between the characteristics of the bioelectrical activity and the EEG.

This paper presents the analysis of the inverse problem for dipolar sources in the brain from measurement of the EEG on the scalp. This analysis uses a simplification in which the original problem is reduced to a problem in a homogeneous region with a null Neumann condition along with a "measurement" (which is obtained from the EEG) on the boundary of the mentioned homogeneous region (Cauchy data). For the analysis and the simplicity of the exposition, the head is modelled for two concentric circles. This allows the forward problem calculation in exact form. Two auxiliary problems are solved in exact form too. For that, the Green function for the Poisson equation with a null Neumann boundary condition in a circular homogeneous region and the circular harmonics are used.

Mathematical Model

We will suppose that the human head, considered as conductive medium, is divided in two disjoint zones as illustrated in Fig. 1. Where Ω = Ω₁ ∪ Ω₂ represents the head, Ω₁ the brain, Ω₂ the rest of the head, σ₁ and σ₂ are the constant conductivities of Ω₁ and Ω₂, S₁ represents the cerebral cortex and S₂ the scalp.

The study of the Inverse Electroencephalogra phic Problem (IEP) can be made through the following boundary value problem [1, 3, 4, 5]:

(1)

(2)

(3)

(4)

(5)

where f is called the source, u_i = _Ωi, i = 1,2, u represents the electrical potential in Ω, and ∂u_i ∂n_j, i = 1,2 denote the normal derivative of u_i on S_j regarding the normal unitary vector n_j, j = 1,2. The boundary conditions (3)-(4) are called the transmission and the condition (5) is obtained when we consider the conductivity of air equal zero. We will call this problem the Electroencephalographic Boundary Problem (EBP).

From the Green formulas the following compatibility condition is deduced

(6)

In the next section, we will study the identification problem of the function f given in (1) using the boundary value problem (1)-(5) and the additional boundary condition:

Simplification of the Inverse EEG problem to one in a homogeneous region

The simplification is achieved through the following steps:

1. We can find the harmonic function u₂ through the Cauchy data u₂ = V and ∂u₂ ∂n₂ = 0 on the boundary S₂. 

2. Making ψ = _S1 we solve the problem to find a harmonic function u in Ω₁ which satisfies the boundary condition ∂u _ ∂n₁ = ψ. We take the solution which is orthogonal to the constants in order to have uniqueness of the solution. The compatibility condition of ψ is given by ∫ _S1ψ(s)ds = 0, which is deduced from the Green formulas.

3. Now, we consider the following inverse source problem: Find the source f that satisfies the problem

using the additional data on S₁:

(8)

where ϕ = _S1.

From this we see that the IEP can be solved through the previous inverse problem. The problem (7-8) is called the Inverse Electroencephalographic Simplified Problem (IESP).

The development of this section is valid for general regions with sufficiently smooth boundaries.

Statement of the inverse problem for dipolar sources

In this paper we are interested in the case where the source is an epileptic focus. In general, the current distributions describing sources of neural activity are quite intricate. A common simplifying assumption is to consider a current dipole as a source. This model, which is known as the equivalent current dipole, has been shown to accurately depict sources not too deep inside the brain [6] and to permit the associated estimation and accuracy analysis to be carried out fairly simply. More complicated shapes may be approximated by multiple dipoles or multipolar expansions. For simplicity, in this work, we consider only a single dipole of current density. In this case, the source f can be represented in the form [1]

(9)

where p represents the dipolar moment (that determines the intensity and orientation of the dipole) and δ(x - a) is the Dirac delta concentrated in a. In this case we have uniqueness of solution for the inverse problem if the measurement V is known on the whole scalp.

We suppose that measurement V is known. With the ideas presented in the previous section, we can recover the parameters of the dipole source, namely, the dipole moment p and its location a, through the minimum squares functional

(10)

where û is the solution of the problem (7) and g is given by (8) and, of course, û is given in terms of the source f and g in terms of V . Some iterative method must be used for the minimization process. In section 6, we apply the function fmincon of MATLAB which uses a Newton type method. 

The solution of the problem (7) when f is given by (9) is [5]

(11)

where G(x,y) is the Green function of the problem (7). The main inconvenience with the Green function technique is due to the difficulty of finding it when the geometry is not simple. 

In the case of general smooth boundaries, in which we can´t use the Green function technique, the ideas presented in this work can be developed using numerical methods (as the boundary element and finite element) since the solution of the problems presented in the last section must be found numerically. This is not considered in this work.

Developing the ideas in the simple case of concentric circles

In order to validate the ideas presented in section 3, the forward problem must be solved. For simplicity of the exposition we consider the case in which the human head is modelled through two concentric circles. We chose this case due to the solution of the forward problem can be calculated in exact form without using numerical methods which is important since we can build synthetic examples for validating the reduction to one problem defined in a homogeneous region. The results of this section can be extended to the case of concentric spheres but in this case we must calculate numerically some integral associated with the Fourier coefficient of the measurement (solution of the forward problem). The case of complex geometries can be solved with the method presented in this work, but it is necessary to use numerical methods for solving the forward and the inverse problems. But even in this case, the proposed method gives advantages because we need fewer calculations. According to the ideas presented in this work the inverse problem is reduced to one in a circular homogeneous region and for this case we will use the Green function technique.

Green’s function for the problem (7)

The Green function G(x,y) defined in Ω₁ × Ω₁ of the problem (7) satisfies the boundary value problem:

(12)

For the case when Ω₁ corresponds to circle of radius R = R₁, the Green function is given by

(13)

where y is the complex conjugate of y [5].

The Green function for the spherical case can be found in [7, 8].

Solution of the forward problem

In order to find the solution of the forward problem, we consider the following auxiliary boundary value problem

(14)

(15)

(16)

(17)

(18)

where g(x) = (x) = _y=a x ∈ S₁. The solution of the problem (14)-(18) is unique in the orthogonal space to the constants [9].

The solution to the problem (1)-(5) is given by

(19)

Now, we calculate the solution of the problem (13)-(17) which can be expressed in the form:

(20)

From (16) and (17) we get

(21)

(22)

where g_kⁱ,i = 1,2 are the Fourier coefficients of g and g(θ) = ∑ _k=1^∞g_k¹ coskθ + g_k² sin kθ.

From (21) and (22) we find

(23)

Now, using the condition (18) we get

(24)

From (23) and (24)

(25)

(26)

(27)

(28)

Substituting (25)-(28) in (21), we find

(29)

(30)

The coefficients (25)-(30) give the solution to the problem (14)-(18) and the measurement is given by V = ∑ _k=1^∞V _k¹ cos(kθ) + V _k² sin(kθ) where

(31)

(32)

Now we calculate the Fourier expansion of the function g = _y=a.

As a first step we must calculate ▿_yG(x,y). After some calculations we get

(33)

Let p = (p₁,p₂). Then

(34)

Taking into account that the solution is orthogonal to the constants, the terms and were neglected. We denote x₁ = R₁ cos(θ), x₂ = R₁ sin(θ), y₁ = r cos(θ_y) and y₂ = r sin(θ_y). After substituting this in (34) we obtain

(35)

Now, we calculate the Fourier coefficients of g.

(36)

(37)

It is necessary to calculate the following integrals

For that we need the following integrals

It can be seen that

(38)

where we use the notation x = R₁e^iθ, y = r₀e^iθ_y. In effect,

and taking into account that dx = ixdθ, then

where now the differential is complex. From this

We define Φ(y) = -. Note that = where x₀ = y. Since Φ is analytical in Ω₁ we have that

Analogously

(39)

Using trigonometric identities we get

(40)

and

(41)

The restriction of the function w₂ to S₂ gives the solution of the forward problem.

Solution of the inverse problem

For solving the inverse problem we suppose that the measurement is given by V = ∑ _k=1^∞V _k¹ cos(kθ) + V _k² sin(kθ). The first step consists of finding the function u₂ (solution of the Cauchy problem in Ω₂ ) from V. After some calculation we find

(42)

The second step consists of finding the harmonic function u in Ω₁, orthogonal to the constants, that satisfies the Neumann condition = ψ = _S1. This function is given by

(43)

From this, the “measurement” g on S₁ is given by

(44)

Note that in (44) the Fourier coefficients of g are in terms of V _k¹ and V _k². On the other hand, these Fourier coefficients are given in terms of the parameters of the dipolar source p and a. The next step consists of determining, from the measurement V, the parameters of the dipole p and a through the least squares functional (10). In the following, we suppose that V _k¹ and V _k² are known. We have to minimize the functional

(45)

where _k¹ and _k² are the approximations of the Fourier coefficients obtained using (44) when the Fourier coefficients of V have errors and N is chosen appropriately such as [10] and [11] (for two and three dimensions, respectively) for the ill-posedness of the problem due to the numerical instability.

When the measurement is given at a finite number of points on the scalp, we can obtain the measurement on the whole scalp using a regularized interpolation method such as presented in [12].

Simple geometries for the analysis of the inverse electroencephalographic problem have been used. For the case of a dipolar source in [13] a method for finding its location is presented which is based on rational approximations in the complex plane. In [14] an algorithm for finding a dipole from Cauchy data is given considering a spherical model for the head.

---

Table 1: Comparison of the exact and approximated parameters for different initial points.

Parameters	r0	θy	p1	p2	Initial point	Euclidian error
Exact	0.7		0.5	0.5
Approximated	0.6667	1.073	0.5197	0.4164	(0.60,0.785,0.60,0.40)	0.0934
Approximated	0.7000	1.594	0.500	0.500	(0.60,0.784,0.58,0.40)	0.0207
Approximated	0.6994	1.0250	0.3239	0.4829	(0.61,0.748,0.58,0.41)	0.0318
Approximated	0.7416	1.0528	0.4405	0.4324	(0.45,0.698,0.65,0.37)	0.0099
Approximated	0.7269	1.1227	0.5719	0.4986	(0.42,0.698,0.68,0.39)	0.0116
Approximated	0.7072	1.0697	0.4028	0.4997	(0.63,2.094,0.60,0.37)	0.0100

---

Numerical examples

In this section we present synthetic examples in order to show the ideas presented in this work. We consider that R₁ = 1.2, R₂ = 2, σ₁ = 3 and σ₂ = 1. We take as the parameters of the dipolar source p = (p₁,p₂) = (0.5,0.5) and (r₀,θ_y) = (polar coordinates). Using (31), (32), (40) and (41) we calculate the corresponding measurement V . In order to simulate the inherent error of the measurement V ^δ due to the equipment, a random error to the Fourier coefficients of V such that _L2(S2) < δ it is included. In this case we have the following approximation to the function g :

(46)

We made programs in MATLAB for determining the parameters of the dipolar source taking as input data the values given in (46) and using the fmincon function which use a Newton type method. In Table 1, the results obtained are presented. For the election of the initial point of the iterative method we can use additional information about the problem (a priori information). We take δ = 0.1 and N = 10.

The ill-posedness of the problem due to the numerical instability must be taken into account for obtaining a stable solution. 

The ideas presented in this work can be used for more realistic geometries of the head along with numerical methods such as finite element method. In [15] numerical methods for dipoles identification are proposed, namely, the finite element method and the boundary element method.

Conclusions

The problem of determining the epileptic focus from EEG on the scalp has been studied through a model that consider the head as a multilayer conductive medium as well as the quasi static approximation of Maxwell equations which address one boundary value problem that allows establishing relationships between epileptic focus and the EEG. The multilayer model can be reduced to one problem in a homogeneous region with a null Neumann condition. From this, the problem of determining the parameters of the dipolar source is analysed in the homogeneous region mentioned above. This is conceptually and numerically more simple than multilayer case. The proposed method is validated numerically through synthetic examples for the case of concentric circles. This simple case was chosen because we were able to obtain the solution of the forward problem in exact form and this fact is used for validating the proposed method. In a similar way, these examples can be developed for the case of multilayer spherical geometry which is commonly used for the study of this problem. Even more, this can be used for more realistic geometries of the head along with numerical methods such as the finite element method. For the case of a finite number of points on the scalp in which the measurement is given, we must apply stable interpolation methods for obtaining the measurement on the whole scalp for the uniqueness of the solution of the inverse problem. However this is not considered in this work.

References

1. Sarvas J. Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem.  Phys. Med. Biol., 1987; 32(1): 11-22.
2. Nuñez PL.  Electric Field of the Brain. Oxford Univ. Press, New York (USA), 1981.
3. Fraguela A, Oliveros J and Morín M. Mathematical Models in EEG inverse, in: Jiménez Pozo MA, Bustamante J, Djorjevich S, Editors, Topics in the Theory of Approximation II (in Spanish). Scientific Texts, Autonomous University of Puebla (México), 2007; 73-95.
4. Grave R, González S and Gómez CM. The biophysical foundations of the localization of encephalogram generators in the brain. The application of a distribution-type model to the localization of epileptic foci (in spanish).  Rev. Neurol., 2004; 39: 748-756.
5. Fraguela A, Morín M. and Oliveros J. Inverse problem statement for locating the parameters of a neural current dipolar source (in spanish).  Communications of Mexican Mathematical Society , 1999; 25: 41-55.
6. Muravchik CH and Nehorai A. EEG/MEG Error bounds for a static dipole source with a realistic head model.  IEEE Transactions on Signal Processing , 2001; 49(3): 470-484.
7. Sobolev SL.  Partial Differential Equations of Mathematical Physics (in spanish). Addison-Wesley Publishing Company, Inc. (Moscú), 1964.
8. Morín M, Oliveros J, Conde J, Fraguela A, Gutierrez EM and Flores E. Simplification of the inverse electroencephalographic problem to one homogeneous region with null Neumann (in Spanish).  Revista Mexicana de Ingeniería Biomédica, 2013; 34(1): 41-51.
9. Fraguela A, Oliveros J and Grebennikov A. Operational statement and analysis of the inverse electroencephalographic problem.  Revista Mexicana de Física, 2001; 47(2):162-174.
10. Fraguela A, Morín M. and Oliveros J. Inverse electroencephalography for volumetric sources.  Mathematics and Computers in Simulation, 2008; 78: 481-492.
11. Oliveros J, Morín M, Conde J and Fraguela A. A regularization strategy for the inverse problem of identification of bioelectrical sources for the case of concentric spheres.  Far East Journal of Applied Mathematics, 2013; 77: 1-20.
12. Bustamante J, Castillo R and Fraguela A. A regularization method for polynomial approximation of functions from their approximate values at nodes.  Journal of Numerical Mathematics, 2009; 17(2): 97-118.
13. Clerc M, Leblond J, Marmorat JP and Papadopoulo T. Source localization using rational approximation on plane sections.  Inverse Problems, 2012; 28(5): 055018 (24pp).
14. Tsitsas N and Martin P. Finding a source inside a sphere.  Inverse Problems, 2012; 28: 015003 (11pp).
15. Pursiainen S, Sorrentino A, Campi C and Piana M. Forward simulation and inverse dipole localization with the lowest order Raviart-Thomas elements for electroencephalography.  Inverse Problems, 2011; 27: 045003 (17pp).