Analysis of the Inverse Electroencephalographic Problem
for Volumetric Dipolar Sources Using a Simplification
J.J. OliverosOliveros ^{*} M.M. MorínCastillo ^{**} F.A. AquinoCamacho ^{*} A. FraguelaCollar ^{*} ^{*} Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de Puebla. ^{**} Facultad de Ciencias de la Electrónica, Benemérita Universidad Autónoma de Puebla. 
Keywords: inverse electroencephalographic problem, volumetric dipolar sources, Green function.

Correspondencia: 
Palabras clave: problema inverso electroencefalográfico, fuentes dipolares volumétricas, función de Green. 
IntroductionDifferent noninvasive 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 noninvasive 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 ModelWe will suppose that the human head, considered as conductive medium, is divided in two disjoint zones as illustrated in Fig. 1. Where Ω = _{1} ∪ Ω_{2} represents the head, Ω_{1} the brain, Ω_{2} the rest of the head, σ_{1} and σ_{2} are the constant conductivities of Ω_{1} and Ω_{2}, S_{1} represents the cerebral cortex and S_{2} the scalp. The study of the Inverse Electroencephalogra phic Problem (IEP) can be made through the following boundary value problem [1, 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
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 regionThe simplification is achieved through the following steps: 1. We can find the harmonic function u_{2} through the Cauchy data u_{2} = V and ∂u_{2}
∂n_{2} = 0 on the boundary
S_{2}. 2. Making ψ = _{S1} we solve the problem to find a harmonic function in Ω_{1} which satisfies the boundary condition ∂ _ ∂n_{1} = ψ. 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_{1}:
where ϕ = _{S1}. From this we see that the IEP can be solved through the previous inverse problem. The problem (78) 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 sourcesIn 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]
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
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]
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 circlesIn 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 Ω_{1} × Ω_{1} of the problem (7) satisfies the boundary value problem:
For the case when Ω_{1} corresponds to circle of radius R = R_{1}, the Green function is given by
where y [5]. The Green function for the spherical case can be found in [7, 8].
Solution of the forward problemIn order to find the solution of the forward problem, we consider the following auxiliary boundary value problem
where g(x) = (x) = _{y=a} x ∈ S_{1}. 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
Now, we calculate the solution of the problem (13)(17) which can be expressed in the form:
where g_{k}^{i},i = 1,2 are the Fourier coefficients of g and g(θ) = ∑
_{k=1}^{∞}g_{k}^{1} coskθ + g_{k}^{2} sin kθ.
Now, using the condition (18) we get
Substituting (25)(28) in (21), we find
The coefficients (25)(30) give the solution to the problem (14)(18) and the measurement is given by V = ∑ _{k=1}^{∞}V _{k}^{1} cos(kθ) + V _{k}^{2} sin(kθ) where
Now we calculate the Fourier expansion of the function g = _{y=a}. As a first step we must calculate ▿_{y}G(x,y). After some calculations we get
Let p = (p_{1},p_{2}). Then
Taking into account that the solution is orthogonal to the constants, the terms and were neglected. We denote x_{1} = R_{1} cos(θ), x_{2} = R_{1} sin(θ), y_{1} = r cos(θ_{y}) and y_{2} = r sin(θ_{y}). After substituting this in (34) we obtain
Now, we calculate the Fourier coefficients of g.
It is necessary to calculate the following integrals
For that we need the following integrals
It can be seen that
where we use the notation x = R_{1}e^{iθ}, y = r_{0}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_{0} = y. Since Φ is analytical in Ω_{1} we have that
Analogously
Using trigonometric identities we get
and
The restriction of the function w_{2} to S_{2} gives the solution of the forward problem.
Solution of the inverse problemFor solving the inverse problem we suppose that the measurement is given by V = ∑ _{k=1}^{∞}V _{k}^{1} cos(kθ) + V _{k}^{2} sin(kθ). The first step consists of finding the function u_{2} (solution of the Cauchy problem in Ω_{2} ) from V. After some calculation we find
The second step consists of finding the harmonic function Ω_{1}, orthogonal to the constants, that satisfies the Neumann condition = ψ = _{S1}. This function is given by in
From this, the “measurement” g on S_{1} is given by
Note that in (44) the Fourier coefficients of g are in terms of V _{k}^{1} and V _{k}^{2}. 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}^{1} and V _{k}^{2} are known. We have to minimize the functional
where _{k}^{1} and _{k}^{2} 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 illposedness 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.
Numerical examplesIn this section we present synthetic examples in order to show the ideas presented in this work. We consider that R_{1} = 1.2, R_{2} = 2, σ_{1} = 3 and σ_{2} = 1. We take as the parameters of the dipolar source p = (p_{1},p_{2}) = (0.5,0.5) and (r_{0},θ_{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 :
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 illposedness 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.
ConclusionsThe 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
