Stable identification of sources located on the cerebral cortex from EEG over the scalp

Arias Cruz, J. A.; Morín Castillo, M. M.; Oliveros Oliveros, J. J.; Gutiérrez Arias, J. E. Moisés; Arias Cruz, J. A.; Morín Castillo, M. M.; Oliveros Oliveros, J. J.; Gutiérrez Arias, J. E. Moisés

doi:10.31349/revmexfis.69.050702

Servicios Personalizados

Revista

Articulo

Indicadores

Citado por SciELO
Accesos

Links relacionados

Similares en SciELO

Otros
Otros

Permalink

Revista mexicana de física

versión impresa ISSN 0035-001X

Rev. mex. fis. vol.69 no.5 México sep./oct. 2023 Epub 28-Oct-2024

https://doi.org/10.31349/revmexfis.69.050702

Gravitation, Mathematical Physics and Field Theory

Stable identification of sources located on the cerebral cortex from EEG over the scalp

J. A. Arias Cruz¹

M. M. Morín Castillo²

J. J. Oliveros Oliveros³
http://orcid.org/0000-0002-1715-3631

J. E. Moisés Gutiérrez Arias⁴

^¹ Facultad de Ciencias de la Electrónica, Benemérita Universidad Autónoma de Puebla, Av. San Claudio y 18 Sur, Colonia Jardines de San Manuel, 72570, Ciudad Universitaria, Puebla Puebla, México. E-mail: angelarias01@gmail.com

^²Facultad de Ciencias de la Electrónica, BUAP, Av. San Claudio y 18 Sur, Colonia Jardines de San Manuel, 72570, Ciudad Universitaria, Puebla Puebla, México. E-mail: maria.morin@correo.buap.mx

^³ Facultad de Ciencias Físico Matemáticas, BUAP Av. San Claudio y 18 Sur, Colonia Jardines de San Manuel, 72570, Ciudad Universitaria, Puebla Puebla, México.

^⁴ Facultad de Ciencias de la Electrónica, BUAP, Av. San Claudio y 18 Sur, Colonia Jardines de San Manuel, 72570, Ciudad Universitaria, Puebla Puebla, México.

Abstract

In this work, we present a stable algorithm for the inverse problem of identifying cortical sources from electroencephalographic measurements on the scalp. This inverse problem is ill-posed due to the numerical instability it presents, i.e., small changes in the measurements can produce large variations in the location of the sources. A boundary value problem is used to find correlations between the sources and measurements. In the case in which the head is modeled using two concentric spheres, we use spherical harmonics to find the solution to the inverse source cortical problem. To handle the numerical instability of this problem, we use the Tikhonov regularization method and a cut-off of the harmonic expansion series. From numerical tests, we found these parameters with which we get good approximations. Finally, we illustrate the algorithm proposed with synthetic examples.

Keywords: Inverse problem; regularization; ill-posed problem; source identification, finite element method

1. Introduction

Electroencephalography is a well-known, non-invasive technique to investigate the brain, based on recording electrodes located on the scalp of its electrical potential (voltage). The EEG is generated by the activity of big conglomerates of neurons that work simultaneously, called bioelectrical sources, which can be located in the cerebral cortex or the brain’s subcortical layers. The electrical activity of neurons is due to the ionic currents. One of the main applications of electroencephalography is to study evoked potentials, for determining zones of bioelectrical activity associated with a specific stimulus, either for studying cognitive processes or for diagnosing and detecting epileptic foci (in both cases, located on the cerebral cortex or subcortical layers). Another application of the EEG is related to determining sleep anomalies. One crucial point is that epileptic foci can be related to edemas, tumors, and calcifications in the brain. In this work, we consider the inverse electroencephalographic problem (EIP), which consists of determining, from electroencephalographic measurements on the scalp, to identify the bioelectrical sources that generate such measurements. In this context, the sources are large clusters of neurons that work in a coupled manner to generate potentials that can be recorded on the scalp employing an electroencephalograph [¹, ²]. The sources can be located in the cerebral cortex or the cerebral volume. The EIP is an ill-posed problem in Hadamard’s sense [³], which is due to:

Different sources can produce the same measurement.
Presents a numerical instability, i.e., small changes in the measurement can produce large variations in the location of the source.

Regularization methods are used to handle this instability, [³] that allow finding stable approximate solutions from the measurement with error. For the study of the EIP, we will first analyze the so-called direct problem, which consists of determining the measurement when the source is known. In the direct problem, unlike the inverse problem, if two sources are close, the respective measurements will also be close.

As we commented above, the bioelectrical sources can be located in the cerebral cortex or the cerebral volume (subcortical). Since the problem is linear regarding the sources, we can separate the study of the inverse problem in cortical o subcortical sources. More precisely, if we have two sources, one cortical and one subcortical, the EEG produced by the two sources is the sum of the EEG produced for each source. The mentioned linearity corresponds with the superposition principle. For this work, we will consider the case of sources on the cerebral cortex and disregard sub-cortical sources’ activity. Moreover, different works study inverse problems in the case of volumetric sources [⁴-⁷]. It was demonstrated in Ref. [⁸] that the EIP has a unique solution. Thus, the ill-posedness of the EIP is associated with numerical instability.

This work presents a stable algorithm to recover cortical sources from EEG measurements when two concentric spheres model the head. In this case, we can use the expansion of the source and the measurement in spherical harmonics. This simplified model has been widely used for studying the inverse source problem since it gives enough information about the brain’s electrical activity [⁹-¹³]. There are numerically comprehensive direct modeling schemes, including realistic geometries and electrical conductivities, even electrically anisotropic tissue. These models have different advantages and disadvantages regarding speed, accuracy, and interpretability of results [¹²]. To handle the numerical instability, we use two regularization parameters, namely,

The Tikhonov regularization parameter.
A truncation of the expansion series.

Both parameters were chosen by numerical tests.

We show the feasibility of the algorithm with numerical examples.

The work is divided as follows: Section 2 shows the boundary value problem that allows the operational formulation of the EIP; Section 3 shows the solution of the forward problem and the Tikhonov regularized solution of the inverse problem, which is used to the cut-off of the expansion series; Section 4 shows numerical examples of the proposed algorithm. Finally, In Sec. 5, the conclusions are shown.

2. Cortical electroencephalographic boundary problem

The analysis of the inverse source cortical problem is carried out by means of the following boundary problem [¹¹]:

Δu1=0 in Ω1, (1)

Δu2=0 in Ω2, (2)

u1=u2 on S1, (3)

σ1∂u1∂n1=σ2∂u2∂n2+g on S1, (4)

∂u2∂n2=0 on S2, (5)

where Ω=Ω1∪S1∪Ω2⊂Rn, with n = 2 or 3, is a region sufficiently smooth. The positive constants σ ₁ and σ ₂ are the conductivities of the regions Ω₁ and Ω₂, respectively (see Fig. 4). Function g represents the source defined on the interface S ₁, while n ₁ and n ₂ are the unit outward normal vectors on the boundary of Ω₁ and Ω₂, respectively. We denote by u the solution of the problem (1)-(5) in Ω and define ui=uΩi, i=1,2. Δ represents to the Laplace operator, which is also denoted by ∇2. The boundary condition (3) corresponds to the continuity of the potential, and condition (4) corresponds to the jump of the normal flux which is equal to the function g, while condition (5) indicates that the conductivity in the exterior region Ω-C is zero (in some applications it corresponds to the conductivity of the air). Ω₁ represents the brain, Ω₂ to the rest of the layers that make up the head (intracranial fluid, skull, scalp), S ₁ represents the surface of the cerebral cortex, S ₂ to the surface of the scalp (see Fig. 2). The function g represents the activity of sources located on the cerebral cortex. We will call the contour problem (1)-(5) the Electroencephalographic Boundary Value Problem and denote it by (EBP).

Figure 1 Schematic representation of the head in two conductive layers. The cerebral cortex is represented by S ₁, which can be considered as an interface between the brain and the rest of the head.

Figure 2 Exact source g(θ,φ)=R13cos(φ)2(φ) cos(2θ).

Figure 3 Approximated source taking 20 terms in Eq. (9).

Figure 4 Ideal measurement A(g) = V obtained solving the forward problem taking 20 terms in Eq. (23).

2.1. Solubility analysis of EBVP

A necessary and sufficient condition for the classical solution of the EBVP is given by [⁸]:

∫S1gds=0. (6)

Definition 1. Given g on S ₁, the Electroencephalographic Forward Problem (EFP) consists of finding the measurement V=uS2, where u is the solution of the EBVP.

Definition 2. Given a function V defined on S ₂, the Electroencephalographic Inverse Problem (EIP) consists of determining a source g defined on S ₁ such that the solution u of the EFP corresponding to g, satisfies that uS2=V.

2.2. The weak solution

2.2.1. Weak solution of the EBVP

We consider the following spaces:

U=g∈L2(S1):∫S1g dS1=0=L2(S1)/R,

V=u∈H1(Ω):∫Ωu dΩ=0=H1(Ω)/R,

W=v∈L2(S2):∫S2v dS2=0=L2(S2)/R.

Definition 3. Given g∈U, a function u∈H1(Ω), is called a weak solution of the problem (1)-(5)) if it satisfies

σ1∫Ω1∇u1.∇vdΩ+σ2∫Ω2∇u2.∇vdΩ=∫S1gv ds, ∀ v∈H1(Ω). (7)

The condition (6) is also necessary and sufficient for the existence of the weak solution of the EBVP [¹¹]. There is a unique weak solution u∈U that satisfies

∥u∥H1(Ω)≤C∥g∥L2(S1), (8)

where the constant C does not depend of g [¹¹].

We want to emphasize that the compatibility condition (6) is necessary and sufficient for the existence and uniqueness of the both classical and weak solutions [⁸].

2.3. Inverse problem. Operational formulation

We can define the following operator T:U→V, defined by T(g)=u, which is continuous due to (8). The composition of the operator T, which is continuous, and the compact trace operator Tr:V→W, defines the compact operator A:U→W, given by A=Tr∘T, i.e.A(g)=uS2, which is injective. Using this operator we will study the problem of identifying sources defined on S ₁.

The inverse operator A ^-1 is not continuous and it leads to numerical instability, which is the cause of the ill-posedness of the problem. To handle the numerical instability, we will use the Tikhonov regularization method. Since Im(A) is dense in 𝒲, we can apply the Tikhonov regularization [³].

The EIP can be expressed in terms of the operator A in the the following form:

Given a functionV∈Wfind a sourceg∈Usuch thatA(g)= V.

Note that the equation A(g)=V does not have solution if V∉Im(A). We can use Tikhonov regularization since Im(A)is dense in 𝒲 [⁸, ¹⁴].

3. EIP Solution

3.1. Solution of the EFP

The boundary value problem (1)-(5) allows to get correlations between EEG and sources. Through this model, we can get a theoretical (exact or ideal) EEG containing the main characteristics to distinguish a real EEG. The boundary value problem (1)-(5) allows to get correlations between EEG and sources. Through this model, we can get a theoretical (exact or ideal) EEG containing the main characteristics to distinguish a real EEG. To achieve this, we solve the so-called forward or direct problem, which consists of determining the exact EEG when the source is known, i.e., when we know the domain and correspondence rule of the function with which the source is represented. To solve the direct problem, as a first step, we find the solution of the EBVP. Then we restrict this solution to the region’s boundary (more precisely, we find the trace of the solution to the region’s boundary).

To find the solution to the problem, it is proposed that the function g∈U is given in the form

g(θ,φ)=∑n=0∞‍∑m=-nn‍gnmYnm(θ,φ), (9)

where the Fourier coefficients gnm of the source g are given by

gnm=∫S1g(θ,φ)Ynm(θ,φ)dS1. (10)

The solution of the EBVP (1)-(5) is sought in the form:

u1=∑n=0∞‍∑m=-nn‍Anm1rnYnm, (11)

and

u2=∑n=0∞‍∑m=-nn‍(Anm2rn+Bnm2r-1-n)Ynm, (12)

where the coefficients Anm and Bnm of the potentials are unknown and must be determined from the boundary conditions of the problem. From the condition u1=u2, on S ₁, it is found that

Anm1R1n=Anm2R1n+Bnm2R1-1-n. (13)

From the condition σ1(∂u1/∂n1)=σ2(∂u2/∂n2)+g on S ₁, it is obtained:

σ1Anm1nR1n-1=σ2Anm2nR1n+Bnm2(n-1)R1-2-n+gnm. (14)

From the boundary condition σ2(∂u2/∂n2)=0, it is found:

σ2Anm2nR2n+Bnm2(n-1)R2-2-n=0. (15)

By subtracting Anm2 from the Eq. (15) we obtain:

Anm2=Bnm2(n+1)R2-2-n/nR2n-1=(n+1/n)Bnm2R2-2n-1. (16)

Multiplying by R ₁ the Eq. (16), we get:

σ1Anm1nR1n=σ2Anm2nR1n+Bnm2(-n-1)R1-1-n+gnmR1, (17)

from where

Anm1=(n+1)R12n+1+nR22n+1nR1n-1Hn gnm. (18)

Multiplying the Eq. (14) by nσ1 yields:

nσ1Anm1R1n=nσ1Anm2R1n+Bnm2R1-1-n. (19)

From (18) and (19), we have the following equation

n(σ1-σ2)Anm1R1n+nσ1+(n+1)σ2×Bnm2R1-1-n=gnmR1. (20)

Substituting Anm2 given by (16), we obtain

Bnm2=R1n+2R22n+1Hn gnm, (21)

where Hn=(n+1)(σ1-σ2)R22n+1+nσ1+(n+1)σ2R2n+2. Substituting Bnm2 on Eq. (16) is obtained

Anm2=R1n+2nHn gnm. (22)

Therefore, the solution of the forward problem for g∈U, is given by

A(g)=∑n=0∞‍∑m=-nn‍(2n+1)R1n+2R2nnEn+FngnmYnm, (23)

where En=n+1σ1-σ2R12n+1, Fn=nσ1+(n+1)σ2R2n+2. Recall that the exact (ideal) measurement V (EEG) is obtained through the relation V(θ,φ)=A(g)(θ,φ).

3.2. Solution of the EIP

From the expression (22) para V=u2S2, we have that the Fourier coefficients of V are given by

Vnm=(2n+1)R1n+2R2nnEn+Fngnm, (24)

from where the coefficients of the source g are given by

gnm=nEn+Fn(2n+1)R1n+2R2nVnm. (25)

Recall that the harmonics Ynm(θ,φ) form an orthogonal basis for the harmonic functions orthogonal to constants, a space arising from the compatibility condition of the function g.

3.3. Solution of the IEP using the Tikhonov regularization method and a cut-off

We consider the equation A(g)=V. To handle the numerical instability of the equation, the Tikhonov regularization method is used, which consists of adding a penalization term to the least square problem J(g)=(1/2)∥A(g)-V∥L2(S2)2, i.e., the method minimizes the Tikhonov functional

Jα(g)=12A(g)-VL2(S2)2+α2∥g∥L2(S1)2, (26)

where α(δ) > 0 and can be chosen using the L-curve method [¹⁵] or the Morozov discrepancy principle [³].

Now, we consider the equation A(g) = V _δ , where ∥V-Vδ∥L2(S2)2, and the functional Jα(g)=(1/2)∥A(g)-Vδ∥L2(S2)2+(α/2)∥g∥L2(S1)2. Let gα(δ)=RαVδ<δ be an approximation of g where Rα=(A*A+αI)-1A*Vδ is the operator associated with the Tikhonov regularization strategy, which has the following properties described in the following theorems [³].

Theorem 1. Rα:W→U is boundedly invertible, and gα(δ)=RαVδ is the unique solution of the normal equation (A*A+αI)gα(δ)=A*Vδ. Furthermore, ∥Rα∥≤1/α.

Theorem 2. Any choice α(δ) such that limδ→0α(δ)=0 and limδ→0δ2 / α(δ)=0 ensures that

supRαVδ-A-1VU : Vδ∈W, V-VδW≤δ→0

when δ → 0. In this case, R _α is called ‘admisible’ and

limδ→0gα(δ)=g in U. (27)

The minimum of the functional (26), with V _δ instead V, is solution of the normal equations [³]:

A*A+αIgα(δ)=A*Vδ, (28)

where A ^* represents the adjoint operator, which is defined through the adjoint boundary value problem

Δw1=0 in Ω1, (29)

Δw2=0 in Ω2, (30)

w1=w2 on S1, (31)

σ1∂w1∂n1=σ2∂w2∂n2 on S1, (32)

σ2∂w2∂n2=ψ on S2, (33)

by the correspondence rule A*(ψ)=σ1w1 with ψ∈W. If

ψ(θ,φ)=∑n=0∞‍∑m=-nn‍ψnmYnm(θ,φ), (34)

it is found that

A*(ψ)=-∑n=1∞‍∑m=-nn‍(2n+1)R2n+2rnnGn×ψnmYnmθ,φ, (35)

where

Gn=n+1σ1-σ2R12n+1+nσ1+n+1σ2R22n+1.

After solving the normal equations, we found that the Fourier coefficients of the regularized source fα(δ) are given by

fnmα=n(2n+1)2n+3GnR2n+2-σ1Kn+αn22n+3Gn2Vnmδ, (36)

where Vnmδ are the Fourier coefficients of V _δ , and Kn=2n+12R12n+3R22n+2. Thus, by (36) we get a stable algorithm to recover the source from the measurement V _δ , which is associated to the EEG.

After numerical tests, we will use 10 terms of the expansion series since we get a good approximation of the regularized sources. Thus, we are using N and α as regularization parameters.

4. Numerical examples

In this section, we present numerical examples to illustrate the algorithm given in Eq. (36). In these examples, we proceed as follows:

The EBVP (1)-(5) is solved for the exact source f.
The forward problem (the trace of the solution of the EBVP) is solved to get the exact or ideal measurement V. We take N = 20 terms in the expansion series.
To get V _δ , we aggregate a random error to the exact measurement V using the rand function of MATLAB to emulate the error in the measurement.
The regularized source fα(δ) is obtained using (36) taking N = 10 and α = 10^-6.
We calculate the relative error given by

RE=f-fα(δ)L2(S1)fL2(S1).

Example 1. We consider the source g(x,y,z)=(x2-y2)z, which in spherical coordinates g(R1,θ,φ)=R13cos(φ)2(φ)cos(2θ). This function satisfies the compatibility condition (6). Figure 2 shows the graphic of this source, and Fig. 3 shows the same graphic obtained using N = 20 terms of the expansion in spherical harmonics. The solution of the forward problem A(g) is shown in Fig. 4, which corresponds to the exact (ideal) measurement taking N = 20 in (23). The error measurement, in Fig. 5, was obtained by adding a random error to each coefficient (using the rand function of MATLAB), taking δ = 0.1. Figure 6 shows the recovered source without regularization and taking N = 90 expansion terms. The relative error between the exact and recovered sources is 9.10. Figure 7 shows the regularized source with α = 10^-6 and N = 10. The relative error is 0.08.

Figure 5 Measurement with error obtained adding a random error to the exact measurement A(g) = V.

Figure 6 Recovered source without regularization (α = 0) and N = 80 terms.

Figure 7 Regularized source obtained taking α = 10^-6 and N = 10 terms in Eq. (36).

Example 2. We consider the source g(x,y)=xz+yz+xy. This function also satisfies the compatibility condition (6). Figure 8 shows the graphic of the exact source. Figure 9 shows the solution A(g) = V of the forward problem. The error measurement is obtained as in the previous example, and it is shown in Fig. 10. Figure 11 shows the recovered source without regularization (α = 0) and N = 80 terms. The relative error is 2.03. Figure 12 shows the regularized source obtained by taking N = 10 and α = 10^-6. The relative error is 0.09.

Figure 8 Exact source g(x, y, z) = (x ² - y ²)z.

Figure 9 Ideal measurement A(g) = V obtained solving the forward problem taking 20 terms in Eq. (23).

Figure 10 Measurement with error obtained adding a random error to the exact measurement A(g) = V.

Figure 11 Recovered source without regularization (α = 0) and N = 80 terms.

Figure 12 Regularized source obtained by taking N = 10 and α = 10^-6.

5. Conclusion

We study the inverse electroencephalographic problem for cortical sources (IEP), which consists of finding the cerebral cortex sources from EEG scalp measurements. This problem is ill-posedness since it presents numerical instability. We present an algorithm to solve the IEP, which uses the Tikhonov regularization and a cut-off (truncation) to handle the numerical instability. We found the regularization parameters from the numerical test to get the stable algorithm. The former regularization parameter is the Tikhonov regularization parameter. The latter is the term N in which we truncate the harmonic expansion of the source. From the numerical tests, we found these parameters, which give good approximations of the sources. In the case of complex geometries for the head, also we can use the Tikhonov functional (26) to recover the source. In this case, it is possible to use an iterative method to find the minimum of the mentioned Tikhonov functional. One of the most used is the conjugate gradient method. However, other methods can be used, like the Quasi-Newton methods. Furthermore, we must solve the direct problems in each iteration using a numerical method, like finite differences or the Finite element method. The parameter associated with the truncation gives an idea to choose the subspace obtained in the discretization of the numerical method used for solving the direct problems. Finally, we illustrate the proposed algorithm with synthetics examples that show its feasibility.

References

1. S. R. Nunez, PL, Electric Field of the Brain, vol. 1 (Oxford University Press: New York, 2006). [ Links ]

2. R. Plonsey, Bioelectrical Phenomena, vol. 1 (Mac Graw Hill, 1969). [ Links ]

3. A. Kirsch, An introduction to the mathematical theory of in verse problems, vol. 120 (Springer, 2011). [ Links ]

4. E. B. A and H. D. T, Some remarks on the problem of source identification from boundary measurements, Inverse Problems 14(1998) 1. [ Links ]

5. M. Morín-Castillo, et al., Simplification of the inverse electroencephalographic problem to one homogeneous region with null Neumann condition, Rev. Mex. Ing. Bio. 34 (2013) 41. [ Links ]

6. A. Amir, Uniqueness of the generators of brain evoked potential maps Trans. Biomed. Eng. 41 (1994) 1. [ Links ]

7. M. Morín-Castillo, et al., Analysis of dipolar sources in the solution of the Electroencephalographic Inverse Problem, Mathematics. 9 (2022) 1. [ Links ]

8. M. Morín-Castillo, et al., Stable identification of sources located on separation interfaces of two different homogeneous media., Adv. Differ. Equ. Control Process. 2019 (2019) 53. [ Links ]

9. N. Tsitsas and P. Martin, Finding a source inside a sphere, Inverse Problems 28 (2012) 1 [ Links ]

10. A. Papargiri, et al., Revisiting an analytical solution for the three-shell spherical human head model in electroencephalography, Partial Differential Equations in Applied Mathematics 4 (2021) 1 [ Links ]

11. A. Fraguela, et al., Operational statement and analysis of the inverse electroencephalographic problem, Rev. Mex. Fis. 47 (2001) 162. [ Links ]

12. S. Næss, et al., Corrected Four-Sphere Head Model for EEG Signals, Frontiers in human neuroscience 11 (2017) 1. [ Links ]

13. F. Vatta, et al., Realistic and Spherical Head Modeling for EEG Forward Problem Solution: A Comparative Cortex-Based Analysis, Comput Intell Neurosci. 2010 (2010) 1. [ Links ]

14. J. J. Conde-Mones, et al., Stable Numerical Identification of Sources in Non-Homogeneous Media, Mathematics. 10 (2022) 1. [ Links ]

15. H. PC and P. D, The use of the L-Curve in the Regularization of Discrete Ill-Posed Problems, BIT Numerical Mathematics 30 (1990) 658 [ Links ]

16. M. Morín-Castillo, et al., Stable Identification of Sources Associated with Epileptic Focus on the Cerebral Cortex, Rev. Mex. Ing. Bio. 40 (2019) 1 [ Links ]

17. J. J. Conde-Mones, et al., Stable Identification of Sources Located on Interface of Nonhomogeneous Media, Mathematics. 9 (2021) 1. [ Links ]

Received: March 07, 2023; Accepted: April 11, 2023

This is an open-access article distributed under the terms of the Creative Commons Attribution License