Analysis of the Inverse Electroencephalographic Problem

for Volumetric Dipolar Sources Using a Simplification

J.J. Oliveros-Oliveros *
M.M. Morín-Castillo **
F.A. Aquino-Camacho *
A. Fraguela-Collar *
* 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.
ABSTRACT

Objective: To analyze the parameter identification problem for volumetric dipolar sources in the brain from measurement of the EEG on the scalp using a simplification which reduces the multilayer conductive medium problem to one homogeneous medium problem with a null Neumann boundary condition.Methodology: The minimum squares technique is used for parameter identification of the dipolar sources. The simple case in which the head is modelled by concentric circles is developed. This case was chosen because we were able to obtain the solution of the forward problem in exact form and for the simplicity of the exposition. Results: The parameter of the dipolar sources can be identified from the EEG on the scalp using the simplification. For the theoretical analysis the results developed for one homogeneous region are used. The numerical implementation is simpler than the multilayer case and the numerical computation requires minor computational cost. Conclusion: The feasibility for solving the parameter identification problem using the simplification is shown. These results can be extended to the case of concentric spheres and complex geometries but the solution cannot be found in exact form.

Keywords: inverse electroencephalographic problem, volumetric dipolar sources, Green function.

Correspondencia:
José Jacobo Oliveros Oliveros
Edificio 111A/201F, 18 Sur y Ave. San Claudio, Col. Jardines de San Manuel, Puebla, Pue. México, CP 72570.
Tel. (01 222) 22 9 55 00 ext. 2119, 2178.
Correo electrónico: oliveros@fcfm.buap.mx
Fecha de recepción:
19 de diciembre de 2013 Fecha de aceptación:
22 de mayo de 2014

RESUMEN

Objetivo: Analizar el problema de identificación de los parámetros para fuentes dipolares volumétricas en el cerebro a partir de la medición del EEG en el cuero cabelludo mediante una simplificación que reduce el problema de un medio conductor de múltiples capas a un problema en un medio homogéneo con una condición de Neumann nula en su frontera.  Metodología: Se utiliza la técnica de mínimos cuadrados para identificar los parámetros de las fuentes dipolares. Se desarrolla el caso simple en el que la cabeza está modelada por círculos concéntricos debido a que la solución del problema directo se puede calcular en forma exacta y por la sencillez de la exposición.  Resultados: Se identifican los parámetros de la fuente dipolar a partir del EEG sobre el cuero cabelludo usando la simplificación. Para el análisis teórico se utilizan los resultados desarrollados para una región homogénea. La implementación numérica es más simple y el cálculo numérico requiere menor costo computacional.  Conclusión: Se muestra la factibilidad para resolver el problema de identificar los parámetros de una fuente dipolar por medio de la simplificación. Los resultados pueden ser extendidos al caso de esferas concéntricas y al de geometrías complejas pero la solución del problema directo no puede hallarse en forma exacta.

Palabras clave: problema inverso electroencefalográfico, fuentes dipolares volumétricas, función de Green.

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 Ω = Ω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, S1 represents the cerebral cortex and S2 the scalp.

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

Δu1 = f   in  Ω1,
(1)

Δu2 = 0   in  Ω2,
(2)

u1 = u2 on  S1,
(3)

  ∂u1     ∂u2
σ1∂n--= σ2∂n--  on  S1,
    1        1
(4)

∂u2
∂n--= 0,  on  S2,
  2
(5)

where f is called the source, ui = u |Ωi, i = 1,2, u represents the electrical potential in Ω, and ∂ui ∂nj, i = 1,2 denote the normal derivative of ui on Sj regarding the normal unitary vector nj, 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

∫
 Ω f(x)dx = 0.
  1
(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:

u2 = V.

PIC
Figure 1. Representation of the head as two homogeneous conductive layers.

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 u2 through the Cauchy data u2 = V and ∂u2 ∂n2 = 0 on the boundary S2

2. Making ψ =      |
σ2∂u2||
σ1∂n1S1 we solve the problem to find a harmonic function u in Ω1 which satisfies the boundary condition u _ ∂n1 = ψ. 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

Δû = f in Ω1,
û ∂n = 0 on S1, (7)

using the additional data on S1:

^u|S1 = g = ϕ - u¯|S1 ,
(8)

where ϕ = u2|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]

f =-1 div (pδ(x - a)),
   σ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

             2
mpi,n a ∥^u(x )- g∥ ,
(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]

       [            ]|
^u(x) =  p-⋅▿  G(x,y) ||  ,
        σ1   y       |y=a
(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 Ω1 × Ω1 of the problem (7) satisfies the boundary value problem:

ΔxG(x,y) = δ(x - y) --1--
πR21 in Ω1,
∂G (x,y)
--------
  ∂nx = 0 onΩ1. (12)

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

           {                        2    2 }
G(x,y) =-1- ln|x- y|+ ln|¯yx - R21|- |x|+-|2y|- ,
        2π                           2R1
(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

Δw1 = 0  in  Ω1,
(14)
Δw2 = 0  in  Ω2,
(15)
w1 = w2 - g on  S1,
(16)
  ∂w1-    ∂w2-
σ1∂n1 = σ2 ∂n1  on  S1,
(17)
∂w2-= 0  on  S ,
∂n2           2
(18)

where g(x) = ^u(x) = [P-           ]||
 σ1 ⋅∇yG (x,y) |y=a x S1. 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

       ({ ^u + w1    x ∈ Ω1
u(x) =
       ( w2        x ∈ Ω2
(19)

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

           ∑∞  1 k        1 k
w1(r,θ)  =     akr coskθ + bkr  sinkθ,
           k=1
           ∑∞ ( 2 k   2 -k)
w2(r,θ)  =      akr + bkr   coskθ
           k(=12 k   2 -k)
        +   ckr + dkr   sinkθ.                               (20)

From (16) and (17) we get

a1Rk  =   a2Rk+ b2R -k- g1,
 k1 1k       k2 1k   k2  1-k   k2
bkR1  =   ckR 1 + dkR 1 - gk,                             (21)

            (           )
σ1a1kRk1 = σ2 a2kRk1 - b2kR-1 k ,
 σ1b1kRk1 = σ2 (c2kRk1 - d2kR-1 k) ,                              (22)

where gki,i = 1,2 are the Fourier coefficients of g and g(θ) = k=1gk1 cos + gk2 sin .

From (21) and (22) we find

         k 2            -k 2     1
(σ1 - σ2)R1ak + (σ1 + σ2)R 1 bk = σ1gk,
 (σ1 - σ2)Rk1c2k + (σ1 + σ2)R -1kd2k = σ1g2k.                           (23)

Now, using the condition (18) we get

a2kRk2 - b2kR-2k = 0,
 c2Rk - d2R -k= 0.                                     (24)
 k  2   k 2

From (23) and (24)

              σ g1Rk
a2k = ---------21k-k-1--------2k,
     (σ1 - σ2)R1 + (σ1 + σ2)R 2
(25)

c2 = ---------σ1g2kRk1---------,
k   (σ1 - σ2)R21k+ (σ1 + σ2)R22k
(26)

2          σ1g1Rk1R22k
bk = (σ---σ-)R2kk+-(σ-+-σ-)R2k,
      1   2  1     1   2  2
(27)

               2 k 2k
d2k =--------σ1gkR-1R2---------.
    (σ1 - σ2)R21k + (σ1 + σ2)R2k2
(28)

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

               (         )
 1  -------σ1g1k-R2k1-+-R22k--------   g1k-
ak = Rk [(σ1 - σ2)R2k+ (σ1 + σ2)R2k] - Rk,
      1          1            2      1
(29)

                (        )
 1         σ1g2k R2k1 + R22k           g2k
bk = Rk[(σ1 --σ2)R2k+-(σ1 +-σ2)R2k] - Rk.
      1          1            2      1
(30)

The coefficients (25)-(30) give the solution to the problem (14)-(18) and the measurement is given by V = k=1V k1 cos() + V k2 sin() where

                  k k
V1k = [--------2σ1R1R2---------]gk1,
     (σ1 - σ2)R2k1 + (σ1 + σ2)R2k2
(31)

 2   [--------2σ1Rk1Rk2--------] 2
Vk =  (σ1 - σ2)R2k1 + (σ1 + σ2)R22k gk.
(32)

Now we calculate the Fourier expansion of the function g = [ p          ]||
 σ1 ⋅▿yG (x,y) |y=a.

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

                        x- y     y
            ∇yG (x,y) = |x--y|2-+ R2-+
((              2)                 1
--x1y1 +-x2y2 --R-1-x1 +-(y1x2---y2x1)x2-,
 (             |x-) y|2             )
 -x1y1 +-x2y2 --R21-x2 --(y1x2 --y2x1)x1
               |x- y|2              .                           (33)

Let p = (p1,p2). Then

           { (              )
g(x)  =   p1  -x1y1 +-x2y2 --R21-x1
         σ1        |x- y|2
          (y x  - y x )x + (x - y )}
     +    -1-2---2-1--22---1---1--
           { (    |x- y|    )
         p2  -x1y1 +-x2y2 --R21-x2
     +   σ2        |x- y|2
                                  }||
     -    (y1x2 --y2x1)x1 +-(x2 --y2) ||                        (34)
                   |x - y|2           y=a

Taking into account that the solution is orthogonal to the constants, the terms y1-
R21 and y2-
R21 were neglected. We denote x1 = R1 cos(θ), x2 = R1 sin(θ), y1 = r cos(θy) and y2 = r sin(θy). After substituting this in (34) we obtain

         p (1- R2 )[R  cosθ- rcosθ ]
g(x)  =  -1------1- --1-------2---y
           ( σ1   )[     |x - y|    ]
      +  p2-1--R21- R1-sin-θ--rsin-θy .                         (35)
             σ2         |x- y|2

Now, we calculate the Fourier coefficients of g.

 1
gk  =  ⟨g(x(),cos(k)θ)⟩
        p1-1--R21-R21∫ 2πcosθ-cos(kθ)
    =        σ1      0    |x- y|2  dθ
          (    2)  2∫ 2π
    +   p2-1--R1-R-1-   sinθcos(k2θ)dθ
       [     σ2     (0    |x- y|        )
    -   (1 - R2)rR   p1cos(θy) + p2sin(θy)
              1   1     σ1         σ2
        ∫ 2π cos(kθ)  ]
    ×       |x-- y|2dθ ,                                         (36)
         0

 2
gk  =  ⟨g(x(),sin(k)θ)⟩ ∫
        p1-1--R21-R21- 2πcosθ-sin(kθ)
    =        σ1      0    |x- y|2  dθ
        p (1- R2)R2 ∫ 2πsinθsin(kθ)
    +   -2-----1---1-   --------2-dθ
       [     σ2     (0    |x- y|        )
    -   (1 - R2)rR1  p1cos(θy) + p2sin(θy)
              1         σ1         σ2
        ∫ 2π sin(kθ)  ]
    ×       |x-- y|2dθ .                                         (37)
         0

It is necessary to calculate the following integrals

     ∫ 2π
I11k =     cos(θ)cos(kθ)dθ,
      0    |x - y|2
 2   ∫ 2πsin-(θ)cos(kθ)
I1k =  0     |x - y|2   dθ,
     ∫ 2πcos(θ)sin(kθ)
I12k =     --------2---dθ,
     ∫0     |x - y|
I2 =   2π sin(θ)sin(kθ)dθ.
 2k   0     |x- y|2

For that we need the following integrals

     ∫ 2πcos(mθ)
Im1 =     ------2dθ,  m = 1,2,3,...
      0  |x- y|

 m   ∫ 2πsin(m-θ)
I2 =  0  |x- y|2dθ,  m = 1,2,3,...

It can be seen that

      2πrm sin(mθ )
Im2 = --m-1---2--y2-  m = 1,2,3,...
     R 1  (R 1 - r0)
(38)

where we use the notation x = R1e, y = r0ey. In effect,

             ( ∫ 2π   m     )
Im2 =  -m1-1Im       --y---2dθ ,
      R1        0  |x- y|

and taking into account that dx = ixdθ, then

               (                  )
 m      --1--      ∫ 2π xm--1--
I2  =   Rm- 1Im   - i 0  |x - y|2ixdθ
         1     (   ∫              )
        --1--             -xm-1--
    =   Rm1- 1Im   - i |x|=R1 |x - y|2dx ,

where now the differential is complex. From this

               (   ∫           m-1      )
Im  =   --1--Im   - i      ---x--------dx
 2      Rm1- 1       |x|=R1 (x- y)(x¯- ¯y)
               (   ∫           m-1       )
    =   -m1--1Im   - i      ---x----x----dx
        R1          |x|=R1 (x- y)(x¯- ¯y)x
          1    (  ∫           xm         )
    =   -m--1Im   i      -------------2-dx
        R1     (   |x|=R1(x - y)(x¯y- R 1) )
          1       ∫           xm
    =   Rm--1Im   i      (y--x)(R2---x¯y)dx  .
         1         |x|=R1        1

We define Φ(y) = ---x2m---
(R1-x¯y). Note that ----xm2----
(x- y)(R1-x¯y) = Φx(-xx)
   0 where x0 = y. Since Φ is analytical in Ω1 we have that

Im2   =  --1--Im (i(2πiΦ(x0)))
        Rm1-1   (         )
        --1--       --ym---    2πrm-sin(m-θy)--
     =  Rm -1Im  2π R21 - r20 =  Rm- 1(R2 - r2).
          1                     1     1   0

Analogously

         m
Im1 = -2πmr-1cos(mθy)-  m = 1,2,3,...
     R 1  (R21 - r20)
(39)

Using trigonometric identities we get

        πrk-1 [r2cos((k + 1)θy)
I11k =   --0k-  -0-(R2--r2)----
         R21        1   0]
    +    R1-cos((k--1)θy) ,
            (R21 - r20)
 2      πrk0-1 [r20sin ((k+ 1)θy)
I1k =   -Rk--  ---(R2---r2)---
          21        1   0]
    -    R1-sin(2(k-21)θy) ,
            (R1[ - r0)
 1      πrk0-1  r20sin-((k+-1)θy)
I2k =    Rk1      (R21 - r20)
         R2 sin((k- 1)θ )]
    +    -1---2---2---y- ,
            (R1[ - r0)
I2  =   πrk0-1  R21cos((k--1)θy)
 2k      Rk1      (R21 - r20)
         r20cos((k + 1)θy)]
    -    ---(R2--r2)---- ,
              1   0
from where
        (    2)    k-1
g1  =   -1--R1-p1πr0---
 k      σ1Rk1-2(R21 - r20)
    ×   [r20cos((k+ 1)θy)+ R21cos((k- 1)θy)]
        (    2)    k-1
    +   -1--R1-p2πr0---
        σ2Rk1-2(R21 - r20)
    ×   [r20sin((k+ 1)θy) +R21sin ((k - 1)θy)]
        [(p cos(θ)   p sin(θ ))
    -     -1----y- + -2----y-
         (   σ1)        σ2    ]
         1---R21-2πrk0+1cos(kθy)
    ×        Rk1-2(R21 - r20)    ,                                (40)

and

        (1- R2)p πrk-1
g2k  =   ---k-12--12-0-2-
        σ[1R1  (R 1 - r0)                ]
    ×   r20sin((k+ 1)θy) +R21sin ((k - 1)θy)
        (1- R2)p πrk-1
    +   ---k-12--22-0-2-
        σ[2R1  (R 1 - r0)                ]
    ×   R21cos((k - 1)θy)- r20cos((k+ 1)θy)
        [(p1cos(θy)   p2sin(θy))
    -     ---σ---- + ---σ----
         (    12)   k+1   2   ]
    ×    1---R1-2πr0--sin(kθy)- .                                (41)
             Rk1-2(R21 - r20)

The restriction of the function w2 to S2 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=1V k1 cos() + V k2 sin(). The first step consists of finding the function u2 (solution of the Cauchy problem in Ω2 ) from V. After some calculation we find

             ∞∑  {   [(   )k   (   )k]
u2(r,θ)  =  1     Vk1   r--  +   R2-   cos(kθ)
           2 k=1       R2       r
               [(   )k  (    )k]      }
        +   V2k   -r-   +  R2-   sin(kθ) .                         (42)
                 R2        r

The second step consists of finding the harmonic function u in Ω1, orthogonal to the constants, that satisfies the Neumann condition -∂¯u
∂n1 = ψ = σ2∂u2||
σ1∂n1|S1. This function is given by

             {        [               ]
       σ2 ∑∞   ( r )k  ( R1)k   (R2 )k
¯u  =   2σ-       R--     R--  -  R--
         1[k=1     1       2   ]}   1
       ×  Vk1cos(kθ)+ Vk2sin(kθ)  .                               (43)

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

         ∞  {(   )  [      ]  (   )  [     ]}
       1∑      R1- k    σ2      R2-k     σ2
g  =   2       R2    1- σ1  +   R1    1+ σ1
        k=[11          2      ]
       × Vk cos(kθ) +V k sin(kθ) .                                   (44)

Note that in (44) the Fourier coefficients of g are in terms of V k1 and V k2. 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 k1 and V k2 are known. We have to minimize the functional

    N
min∑  (g1 - ˜g1)2 +(g2 -g˜2)2,
p,ak=1  k   k      k   k
(45)

where ˜g k1 and ˜g k2 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].


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







Parameters r0 θy p1 p2 Initial point Euclidian error







Exact 0.7 π
3 0.5 0.5
Approximated0.6667 1.073 0.51970.4164(0.60,0.785,0.60,0.40) 0.0934
Approximated0.7000 1.594 0.500 0.500 (0.60,0.784,0.58,0.40) 0.0207
Approximated0.69941.02500.32390.4829(0.61,0.748,0.58,0.41) 0.0318
Approximated0.74161.05280.44050.4324(0.45,0.698,0.65,0.37) 0.0099
Approximated0.72691.12270.57190.4986(0.42,0.698,0.68,0.39) 0.0116
Approximated0.70721.06970.40280.4997(0.63,2.094,0.60,0.37) 0.0100








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 examples

In this section we present synthetic examples in order to show the ideas presented in this work. We consider that R1 = 1.2, R2 = 2, σ1 = 3 and σ2 = 1. We take as the parameters of the dipolar source p = (p1,p2) = (0.5,0.5) and (r0y) = (0.7,π/  )
       3 (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 ∥∥      ∥∥
∥V - V δ∥L2(S2) < δ it is included. In this case we have the following approximation to the function g :

          ∞ { (   )k[      ]   (   )k[      ]}
gδ  =   1∑     R1-   1 - σ2  +  R2-   1 + σ2
        2k=1   R2        σ1     R1        σ1
         [  1,δ          2,δ      ]
       ×  Vk  cos(kθ)+ Vk  sin (kθ) .                                 (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).