Revista Mexicana de Ingeniería Biomédica

 

Simplificación del Problema Inverso

Electroencefalográfico a una Sola Región

Homogénea con Condición de Neumann Nula

INTRODUCCIÓN

En diversos campos de la investigación, se presentan situaciones en las cuales es necesario conocer las causas que producen cierto fenómeno a través de la información parcial que se obtiene del mismo [1]. Este tipo de problemas son llamados de identificación y son ampliamente estudiados en muchos campos de la investigación, entre otros la medicina, donde hay un gran interés en el problema de identificación de fuentes bioeléctricas cerebrales, a partir de los datos obtenidos por medio de un electroencefalograma (EEG) ya que permiten detectar posibles anomalías (daños, mal funcionamiento, etc.) lo que se ha hecho tradicionalmente a través de diferentes técnicas de diagnóstico.

Para los problemas de identificación se usan modelos matemáticos, con los cuales se desarrollan técnicas no invasivas. Entre estas se encuentran la tomografía por emisión de positrones, la resonancia magnética nuclear, la electroencefalografía, la cual es de particular de interés para nosotros. Por medio de esta técnica se registran los potenciales en un electroencefalograma (EEG); estos potenciales provienen de la actividad eléctrica de los tejidos excitables, y se captan midiendo la diferencia de potencial existente entre un electrodo explorador y otro de referencia. A las fuentes que son generadas por la actividad electroquímica de estos órganos se les conoce como fuentes bioeléctricas y se considera que están compuestos por grandes conglomerados de neuronas que actúan simultáneamente. Entre las ventajas de la técnica del EEG se encuentran que la información que proporciona se captura en tiempo real, de manera simple, es no invasiva además de económica. En algunos casos podemos considerar que los generadores están concentrados en una región del cerebro y que pueden representarse por funciones de cuadrado integrable definidas sobre esa región. En el caso particular en que se tiene una fuente dipolar (con la cual se representan a los focos epilépticos), es necesario enfocar el análisis de este problema a través de las distribuciones o funciones generalizadas. Para el estudio de este problema se consideran el problema directo e inverso electroencefalográfico. El primero de ellos consiste en hallar el EEG cuando conocemos a la fuente. El segundo consiste en determinar la fuente a partir del EEG.

En este trabajo se presenta una simplificación del PIE, la cual se basa, como primer paso, en separar el problema en dos sub problemas. El primero de ellos puede ser calculado a partir del EEG en el cuero cabelludo. Para el segundo de ellos, se divide en otros dos, uno de los cuales puede calcularse también (indirectamente) de la medición. De esta manera, el problema inverso se reduce a un problema de identificación de fuentes en un medio homogéneo con una condición de Neumann nula en la frontera de dicho medio, con lo que se pueden usar los resultados que se tienen para este último problema el cual ha sido utilizado para realizar experimentación y para hallar resultados teóricos. Para validar el método se presenta la solución del problema inverso para el caso de fuentes distribuidas en el volumen cerebral. Para fuentes dipolares (focos epilépticos) se halla la solución para una región anular homogénea usando una función de Green. En el caso de una esfera la función de Green puede calcularse en forma explícita y el potencial producido en un punto se obtiene realizando un producto punto de vectores en el espacio por lo que el problema inverso se reduce a un problema de mínimos cuadrados no lineales en un medio homogéneo. En el desarrollo de este trabajo no se considera el caso de errores en la medición.

Debido a que la epilepsia ataca a una parte importante de la población mundial [2] se han dedicado esfuerzos para su estudio los cuales pueden dividirse, en desarrollo de algoritmos y en desarrollo tecnológico para obtener sistemas de visualización de la zona afectada. Este trabajo se enmarca dentro del primero de ellos.

Modelo Matemático

El modelo matemático que se presenta en este trabajo ha sido ampliamente utilizado para el estudio del problema de identificación en ([1], [3], [4], [5]); en él la cabeza humana ha sido modelada por medio de capas conductoras con conductividad constante y diferente en cada capa. La actividad eléctrica del cerebro es registrada en el cuero cabelludo por medio del EEG. Se considera para la modelación que el EEG es producido por grandes conglomerados de neuronas que se activan simultáneamente. A estos conglomerados se les conoce como generadores o fuentes bioeléctricas ([4], [6]). Adicionalmente, se supone que las corrientes que pueden producirse en la región Ω se deben únicamente a la actividad eléctrica del cerebro y pueden ser de dos tipos: óhmicas e impresas. Las primeras se deben al movimiento de cargas iónicas a través del fluido extracelular en el cerebro y las segundas a las corrientes de difusión a través de las membranas neuronales ([6]) las cuales se denotan por J^p y son las de interés en el problema de identificación, ya que el soporte de estas puede darnos información sobre la ubicación espacial de la zona afectada. Tomando en cuenta lo anterior, se demuestra que el estudio del Problema Inverso Electroencefalográfico (PIE) puede ser realizado a través del siguiente problema de valores en la frontera ([1], [4], [5], [7]):

(1)

(2)

(3)

(4)

(5)

donde Ω = Ω₁ ∪ Ω₂ representa a la cabeza, Ω₁ el cerebro, Ω₂ el resto de las capas que componen la cabeza (líquido intracraneal, cráneo, cuero cabelludo), σ₁ y σ₂ son las conductividades de Ω₁ y Ω₂ las cuales se suponen constantes (ver figura 1), f = divJ^p∕σ₁ es llamada la fuente, u_i = _Ωii = 1,2 y u representa al potencial eléctrico en Ω. El símbolo Δ representa al operador laplaciano, que también se simboliza como ∇². Obsérvese que en este caso podremos recuperar sólo una parte de la fuente bioeléctrica J^p. Las condiciones de frontera (3)-(4) son llamadas de transmisión y la condición de frontera (5), se obtiene al considerar que la conductividad de Ω^C es cero (la conductividad del aire). De las fórmulas de Green se deduce la siguiente condición de compatibilidad

(6)

En el caso en que sólo se consideren fuentes corticales, su presencia se refleja sobre la condición de frontera asociada con la igualdad de flujos de corriente. Si denotamos por j^p a la densidad de corriente cortical entonces dicha condición de frontera toma la forma:

Se denomina a j^p ⋅n₁ fuente cortical y satisface una condición de compatibilidad similar a la (6). Sin embargo, si consideramos simultáneamente tanto fuentes corticales como volumétricas esta condición se convierte en ∫ _Ω1f(x)dx + ∫ _∂Ω1j^p ⋅n₁dx = 0. El caso de fuentes corticales ha sido estudiado en [8] y se busca la solución como la suma de un potencial de capa doble definido S₁ más uno de capa simple definido sobre S₂. En este caso la densidad dipolar definida sobre S₁ puede representar la actividad de neuronas piramidales y, por lo tanto, brindar información sobre la zona activa de la corteza cerebral. Sin embargo, en ese trabajo sólo se consideraron fuentes que pueden representarse por funciones de cuadrado integrable y no por funciones generalizadas con las cuales se representan las fuentes dipolares. Este caso debe estudiarse matemáticamente con detalle para asegurar la existencia de soluciones de este problema además de interpretarse apropiadamente el planteamiento matemático para que dichas soluciones tengan un sentido en el problema real. La discusión sobre este tema no se aborda en este trabajo por lo que no se considera la presencia de fuentes corticales.

Al problema (1)-(5) se llama Problema de Contorno Electroencefalográfico (PCE).

Figura 1: Representación de la cabeza como dos medios conductores homogéneos acoplados.

Simplificación del problema a una región

Problemas que permiten la simplificación

Para resolver el Problema Inverso Electro-encefalográfico, primero se desacopla el problema (1)-(5) en los dos problemas siguientes:

Primer problema. Resolver el problema de Cauchy en la región anular Ω₂ dado V sobre S₂:

(7)

Segundo problema. Después de hallar u₂ en el primer problema por medio de los datos de Cauchy sobre S₂ usando las condiciones de transmisión (3)-(4) el PIE se reduce a hallar f a través del problema:

(8)

donde φ = u₂|_S1 y ψ = _S1.

La condición de transmisión (4) correspondiente a la igualdad de los flujos de corriente, lleva a que

(9)

Para reducir a una región, consideremos el

Tercer problema. Para el estudio del problema (8) tomamos los siguientes dos problemas:

(10)

(11)

Notemos que u₁ = û + u si elegimos f en (11) a través del dato adicional:

(12)

La solución del problema (10) se toma ortogonal a las constantes, con lo cual es única. Debido a que ψ es conocida a partir de los datos de Cauchy, el problema (10) puede separarse del resto y sólo considerarse para el estudio del problema de identificación, al problema (11)-(12), el cual corresponde al caso del PIE para una sola región con una condición de contorno de Neumann nula y que llamaremos Problema Inverso Electroencefalográfico Simplificado (PIES). Nótese que la condición de compatibilidad sobre ψ para el problema (10) es la dada por (9). Llamaremos al problema (11), Problema de Contorno Electroencefalográfico Simplificado (PCES).

El desarrollo de esta sección es válido para regiones generales con frontera suficientemente suave que incluyen el caso de las geometrías que pueden aproximar con mayor precisión a la geometría de la cabeza de un paciente específico. En la sección siguiente se desarrollan los tres problemas anteriores para el caso en el que la cabeza se modela por medio de esferas concéntricas, a fin de ilustrar cómo funciona la propuesta realizada en esta sección.

Ejemplificación de la técnica propuesta usando el caso de esferas concéntricas

Es conocido [9] que para fuentes bioeléctricas volumétricas representadas por funciones de cuadrado integrable, la solución del PIE es única si la fuente es una función armónica y satisface la condición (6).

Para ejemplificar la metodología sobre la reducción a una sola región, consideremos el caso de fuentes f representadas por funciones armónicas que satisfacen (6) (con lo que, como se dijo anteriormente, garantizamos unicidad de la solución del problema de identificación) en el caso en que la cabeza se modela por dos esferas concéntricas lo cual es elegido por claridad de la exposición. En este caso se pueden realizar todos los cálculos de manera analítica sin recurrir a los métodos numéricos. Sin embargo, los resultados se pueden generalizar al caso en el que las capas no son necesariamente esféricas para lo cual sería necesario utilizar métodos numéricos para validar los resultados que se presenten.

Primer problema. Resolver el Problema de Cauchy en la región anular Ω₂ compuesta por esferas concéntricas, dada la medición sin error V :

(13)

Sean Y _nm los armónicos esféricos y V _nm los coeficientes de Fourier de V en la base que ellos forman. En este caso se tiene que V está dado por

(14)

de donde la correspondiente solución del problema (13) es

u2(r,θ,ϕ) = ∑ n=0∞∑ m=-nn
V nmY nm(θ,ϕ).		(15)

Debe observarse que para que la serie dada en (15) converja, los coeficientes de Fourier deben converger muy rápidamente debido a la presencia del término ⁿ⁺¹.

El potencial u₂ sobre S₁ es

φ = ∑ n=0∞∑ m=-nn
V nmY nm(θ,ϕ).		(16)

Usando la expresión (15) puede ser calculado sobre S₁:

ψ = S 1 = r=R 1 = ∑ n=0∞∑ m=-nn
V nmY nm(θ,ϕ).		(17)

Segundo problema. Resolver el problema (10):

Como u es armónica esta se busca en la forma u(r,θ,ϕ) = ∑ _n=0^∞∑ _m=-nⁿA_nmrⁿY _nm(θ,ϕ). Se halla que

Utilizando la condición _S1 = _S1 = ψ se tiene que

Por lo tanto, la solución del problema (10) viene dado por

u(r,θ,ϕ) = ∑ n=1∞∑ m=-nn
V nmrnY nm(θ,ϕ)

de aquí que

S1 =	∑ n=0∞∑ m=-nn
	V nmY nm(θ,ϕ).	(18)

Tercer problema. En este paso se resuelve el PIES.

Sea = ∑ _n=0^∞∑ _m=-nⁿ _nmY _nm(θ,ϕ), donde = φ -_S1. De las ecuaciones (16) y (18) se tiene que

de donde los coeficientes de Fourier de son:

(19)

Ahora, el problema inverso de identificación de la fuente f se puede estudiar en la región homogénea Ω₁ y consiste en hallar f a partir del problema:

con la información adicional û = sobre S₁. En este caso la solución del PIES está dada por

(20)

Figura 2. Comparación de las superficies de nivel de las fuentes exacta y recuperada para r = 1.

la cual está en función de los coeficientes de Fourier _nm de dados por (19) (que a su vez están dados en términos de los coeficientes de Fourier de la medición V ). En términos prácticos la serie (20) debe truncarse para poder graficarla. Más aún, para el problema directo las series deben truncarse en el entendido que entre más términos se tomen de la misma, mejor será la aproximación. En el problema inverso el truncamiento debe hacerse en función del error en los datos. Como en este caso no estamos considerando error en las mediciones no se ha analizado este hecho.

En el siguiente ejemplo sintético, la serie (20) asociada a la fuente sólo tiene un número finito de términos.

Ejemplo. Consideremos el caso en que R₁ = 1,R₂ = 1.2,σ₁ = 3yσ₂ = 1 y

En coordenadas cartesianas f(x,y,z) = 2z. Esta fuente produce la medición

(21)

Los coeficientes de Fourier de la medición V son

V 5,4	=
	= V 5,-4.

Usando el método propuesto calculamos :

donde sus coeficientes de Fourier son

Sustituyendo estos coeficientes en la ecuación (20) se halla que la fuente correspondiente a la medición dada por (21) está dada por:

f(r,θ,ϕ) =	-5
	+ .

Las representaciones gráficas de las fuentes exacta y recuperada, están dadas por la figura 2.

Como no se consideran errores en la medición, la fuente exacta y la recuperada prácticamente coinciden. De hecho, en la tabla I se muestran los valores para los diferentes tipos

de error. Mediante EA, ER y EMC denotamos el error absoluto, el error relativo y el error en media cuadrática entre la fuente exacta y la recuperada, respectivamente.

Los valores presentados en la Tabla I, muestran que el método recupera a la fuente exacta considerando aritmética exacta. Para el caso de errores en la medición, debido a que el problema de Cauchy (8) es mal planteado [9], la fuente obtenida por el método puede variar sustancialmente de la fuente exacta y más aún, pequeños cambios en la geometría también pueden provocar también cambios sustanciales en la solución del problema de identificar a las fuentes. Para ver esto, basta considerar que en vez de tener la medición sobre el cuero cabelludo, la tenemos en una curva paralela muy cercana al mismo (hacia el interior de la cabeza). Para tener el valor de la medición en dicha superficie paralela, debemos resolver el problema de Cauchy para la ecuación de Laplace en la región anular formada por estas dos superficies, y por lo dicho arriba sobre su inestabilidad, tendremos cambios importantes en el valor del potencial en la mencionada superficie paralela, y por ende, en la localización de la fuente. Por ello, en este caso es necesario usar métodos de regularización como el de Tijonov que tomen en cuenta lo mencionado arriba.

Planteamiento del problema inverso para fuentes dipolares

En esta sección plantearemos el problema de determinar los parámetros de una fuente dipolar usando la simplificación presentada en este trabajo. Para ello, dada la medición V sobre el cuero cabelludo, calcularemos dada por (12) siguiendo los pasos descritos en las secciones anteriores. Estudiaremos ahora el problema (11).

Solución del PCES

En lo que sigue consideraremos el problema (11), al que hemos denotado por PCES, pero por simplicidad utilizaremos u en vez de û y Ω en vez de Ω₁.

Como primer paso para el estudio del problema de identificación se debe estudiar la solubilidad del PCES. En este trabajo estamos interesados en el caso en que la fuente corresponde a un foco epiléptico y la representación matemática de este tipo de fuentes es por medio de las funciones generalizadas o distribuciones. Más precisamente, un foco epiléptico concentrado en el punto a puede representarse en la forma ([4]):

(22)

donde p representa al momento dipolar y δ(x - a) es la función delta de Dirac concentrada en a. La idea básica para hallar la solución en este caso es la siguiente: ya que δ(x - a) es el límite de funciones suaves concentradas alrededor de a, para cada una de dichas funciones suaves, se busca la solución clásica en términos de una función de Green y se toma el límite de dichas soluciones clásicas para hallar la solución que corresponde a f = divJ^p∕σ₁ cuando J^p está dado por (22). Para ello comencemos con la siguiente

Notemos que una condición necesaria de existencia de la solución clásica del PCES es que:

(23)

la cual se obtiene de las fórmulas de Green. Además la solución clásica del PCE es única salvo constantes y puede expresarse en términos de la función de Green que se define a continuación.

Consideremos los siguientes espacios: C₀(Ω) es el conjunto de funciones continuas que se anulan en una vecindad de la frontera de Ω y C₀^⊥ = .

Para cada f_n ∈ C₀^⊥(Ω), a partir de (22), la solución del PCES, con f sustituido por f_n, está dada por:

(27)

Si además f_n ∈ C₀^⊥(Ω) converge a f = divJ^p∕σ₁, con J^p dada por (19), se tiene que lim _n→∞u_n existe y no depende de la sucesión que converja a f. A tal límite se le llama solución del PCES ([7]).

Pasemos ahora al caso en que la fuente bioeléctrica está dada por

(28)

La solución del PCES cuando J^p = pδ(x - a) viene dada por:

(29)

Este resultado puede consultarse en [7]. De esta forma el potencial producido por una fuente dipolar, se expresa en función del momento dipolar y el gradiente de la función de Green.

Para la solución de este problema pueden usarse otras técnicas como la presentada en [10], en donde se propone un método libre de mallas que es computacionalmente eficiente. Sin embargo, para el caso de fuentes dipolares la expresión (29) para x ∈ S es muy sencilla, representa a la solución del problema directo, es decir, al potencial teórico medido en la frontera de Ω. Claramente es necesario conocer a la función de Green en (29) para poder calcular el potencial teórico. Esto se realiza en la siguiente sección para el caso en que la región Ω es una esfera.

Caso de una región esférica

Función de Green

En esta sección consideramos el caso en que Ω corresponde a una esfera de radio uno en la cual la función de Green está dada por ([11]):

(30)

donde r = |x - y|, g(x,y) = αR² + g₁(y,x),paraα = -, R² = ||y||², m(Ω) es el volumen de Ω y g₁ una función armónica que satisface el problema

(31)

(32)

La elección de α = -1∕6m(Ω) garantiza la existencia y la unicidad (salvo constantes) de la función g₁ y, por lo tanto, de la función de Green G.

Cuando x está sobre el eje z, es decir, cuando x = (0,0,z₀) se tiene que

(33)

en donde

z0′ = , r1 = |y - x′|, x′ = ,
w = log e(z0′- z + r1),C = - -log ez0.

La constante C se elige de la condición:

(34)

con la cual garantizamos que la función de Green es única.

Nótese que cuando z₀ = 1 se tiene que el punto x corresponde al polo norte de la esfera. Tendremos entonces el potencial en ese punto. Para el caso en que x no coincida con el polo norte, se realiza un cambio de coordenadas para colocar el eje z en la dirección del punto x y que éste coincida con el polo norte (para lo cual se usan matrices de rotación). Así tendremos el potencial teórico en cada punto de la superficie.

Modelos que consideran una geometría real de la cabeza [10], han sido empleados para el estudio de este problema cuando se tiene una fuente dipolar, los cuales podrían dar información sobre la posible ubicación de la fuente. Sin embargo, de acuerdo con [12] un modelo muy detallado de la cabeza humana podría no ser necesario para el estudio tanto del problema directo como del inverso.

Planteamiento del problema inverso

A partir de los resultados presentados en este trabajo, el problema de determinar los parámetros de una fuente dipolar, a saber, el momento dipolar p y la posición a, puede estudiarse a partir del funcional de mínimos cuadrados

donde x pertenece a la superficie de la esfera interior que representa al cerebro, u está dada por (29) y la función de Green por (33). Se debe proponer un algoritmo estable el cual puede basarse en la técnica de mínimos cuadrados no lineales planteado en este trabajo junto con algún método de regularización para tratar la inestabilidad del problema.

Conclusiones

En este trabajo por medio de una simplificación se redujo el PIE al análisis de un problema de Poisson (definido en una sola región homogénea con conductividad constante) con una condición de contorno de Neumann nula junto con el dato que se determina de la medición sobre el cuero cabelludo (ver sección 3). Este caso ha sido estudiado como un primer paso para el análisis del PIE. En este trabajo se demuestra que este problema tiene una relación estrecha con el problema de varias capas conductoras. Esto se deriva del hecho que el flujo de corriente es nulo en la superficie del cuero cabelludo debido a que el aire se considera un aislante. En aquellos casos en que esta condición física se cumpla, pueden aplicarse los resultados presentados aquí.

Se validó, usando la técnica de los armónicos esféricos, la simplificación en el caso en que la cabeza se modela por esferas concéntricas y la fuente es una función armónica. Para el caso de fuentes dipolares, se planteó el problema de identificar los parámetros de una fuente dipolar usando la simplificación. Para el caso en que las mediciones tengan errores, se debe proponer un algoritmo regularizado ya este problema es mal planteado debido a que pequeños errores en la medición, pueden generar grandes cambios en la localización de la fuente.

Referencias

1. Fraguela A, Morin M. Oliveros J. “Modelos Matemáticos en Electroencefalografía Inversa” En: Jiménez Pozo M, Bustamante González J, editores, Tópicos de la Teoría de Aproximación II, Capítulo 4. Textos Científicos, Benemérita Universidad Autónoma de Puebla, 2007; 73-95.
2. Fuertes de Gilbert Rivera B., López Gutiérrez R., Gil Gregorio P. “EPILEPSIA” En: Tratado de geriatría para residentes. Editado por Sociedad Española de Geriatría y Gerontología (SEGG), Príncipe de Vergara, 57-59, 1.ž B, 28006 (Madrid), 2006: 519-530.
3. Heller L. “Return Current in Encephalography. Variational Principles”. Biophysical Journal, 1990; 5:. 601-607.
4. Sarvas J. “Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem”. Phys. Med. Biol., 1987, 32(1): 11-22.
5. Grave R, González S, Gómez C. M. “Bases biofísicas de la localización de los generadores cerebrales del electroencefalograma. Aplicación de un modelo tipo distribuido a la localización de focos epilépticos”. Revista de Neurología, 2004; 39: 748-756.
6. Nuñez P.L. “Electric Field of the Brain". New York, Oxford Univ. Press, 1981.
7. Fraguela A, Morín M, Oliveros J. “Planteamiento del Problema Inverso de Localización de los Parámetros de una Fuente de Corriente Neuronal en forma de Dipolo”. Comunicaciones de la Sociedad Matemática Mexicana, 1999; 25: 41-55.
8. Fraguela A, Oliveros J, Grebennikov A. “Planteamiento operacional y análisis del problema inverso electroencefalográfico”. Revista Mexicana de Física, 2001; 47 (2): 162-174.
9. Fraguela A, Morín M, Oliveros J. “Inverse electroencephalography for volumetric sources”. Mathematics and Computers in Simulation, 2008; 78: 481-492.
10. Von Ellenrieder N, Muravchik C. H, Nehorai F. A. “A Meshless Method for Solving the EEG Forward Problem”. IEEE Transactions on Biomedical Engineering, 2005; 52(2): 249-257.
11. Sobolev S. L. “Partial Differential Equations of Mathematical Physics”. Addison-Wesley Publishing Company (Inc. Moscú), 1964.
12. Von Ellenrieder N, Muravchik C. H, Nehorai F. A. “Effects of Geometric Head Model Perturbations on the EEG Forward and Inverse Problems”. IEEE Transactions on Biomedical Engineering, 2006; 53(3): 421-429.