JORNADAS CHILENAS DE INGENIERIA BIOMEDICA 2007

Detección Automática de Husos Sigma en Registros Polisomnográficos de Niños Leonardo Causa1, Claudio M. Held1, Javier Causa1, Carlos Benavides1, Pablo Estévez1, Claudio Pérez1, Marcelo Garrido2, Cecilia Algarín2, Patricio Peirano2. 1 Laboratorio de Ingeniería Biomédica, Departamento Ingeniería Eléctrica, Universidad de Chile 2 Laboratorio de Sueño, Instituto de Nutrición y Tecnología de los Alimentos, Universidad de Chile [email protected]

Resumen—Presentamos un sistema de detección automática de Husos Sigma (HS) en registros polisomnográficos de niños. Se utiliza Descomposición Modal Empírica (EMD) y Transformada de Hilbert-Huang (THH) para determinar la amplitud y frecuencia instantánea de las señales de las derivaciones anteriores (F4-C4 y F3-C3) del EEG. Luego se aplican criterios de experto para generar candidatos a HS, y se consideran aspectos morfológicos y frecuenciales para determinar los HS. El algoritmo fue entrenado sobre un conjunto de trozos de registros de EEG en distintas estados y etapas. El conjunto de prueba consta de 4 registros completos de sueño, equivalentes aproximadamente a 34,2 horas de sueño nocturno. Los resultados muestran un 90,6% de concordancia con la inspección visual del experto y un 74,1% de precisión en la detección. Palabras Claves—Husos Sigma, EEG, Descomposición Modal Empírica (EMD), Transformada de Hilbert-Huang, Detección automática de patrones de sueño I.

E

INTRODUCCION

L estudio del sueño y de los patrones característicos asociados es usado con diferentes propósitos. En especial, el estudio de sueño en niños se relaciona con la evaluación del desarrollo neurofuncional del Sistema Nervioso Central (SNC), deficiencias nutricionales, entre otros. El polisomnograma es un registro de distintas señales fisiológicas adquiridas durante el sueño, tales como el electroencefalograma (EEG), electrooculograma (EOG), electromiograma (EMG), movimientos corporales (MC), electrocardiograma (ECG), oximetría, etc. Se distinguen dos estados de sueño, Sueño Quieto (SQ) y Sueño Paradójico (SP); el SQ se divide en cuatro etapas (SQ14). La evolución del sueño a través de los distintos estados y etapas puede visualizarse a través del hipnograma. La detección visual de patrones de sueño requiere de gran dedicación en tiempo y esfuerzo por parte del experto, además de ser un proceso subjetivo y que presenta variabilidad entre distintos expertos ([1],[2]). La detección automática de estados y patrones de sueño es una poderosa herramienta para la disminución del tiempo y variabilidad en el análisis. Los Husos Sigma (HS) son patrones de sueño que se presentan en forma de ráfagas de actividad altamente rítmica en el rango de frecuencias entre 10-15 Hz que aparecen predominantemente en las derivaciones frontales del EEG. Son de apariencia fusiforme y su duración es mayor a 0,5 s.

Normalmente se encuentran montados sobre ondas más lentas de mayor amplitud. La figura 1 muestra ejemplos de HS en un registro EEG. Los HS se asocian a la determinación del inicio del SQ2 [3], y también a procesos de consolidación de la memoria, principalmente de tipo declarativa. Gais et al. [4] mostraron que la densidad de HS durante la noche aumenta significativamente después de la realización de una tarea de aprendizaje, en comparación con tareas de control. Estos resultados mostrarían que los HS serían sensibles a la experiencia de aprendizaje.

Fig. 1. Ejemplos de 3 Husos Sigma marcados por el experto en derivación F4-C4 (C1) de un registro EEG en etapa SQ2. Distintos grupos han realizado trabajos para la detección automatizada de HS. Huupponen et al. [5] muestran una comparación entre 4 métodos distintos. Los mejores resultados presentan un 70% de detección y un 98,6% de precisión sobre trozos de registros previamente identificados como SQ2, utilizando un método que combina Transformada Rápida de Fourier (FFT) y umbrales de amplitud. Los mismos autores en [6] describen un módulo para establecer un umbral adaptivo para la detección de HS, reportando como mejor resultado un 84,9% de detección en 4 registros completos. Gorur et al. [7] utilizan Short Fourier Transform (SFT) y luego redes neuronales (MLP) y máquinas de soporte vectorial (SVM), obteniendo un 88,7% de detección con MLP y un 95,4% con SVM, en trozos de registro de SQ2. En nuestros trabajos anteriores ([8],[9],[10]) se han utilizado diversas características basadas en criterios de experto para la detección de HS. En [8] se obtuvo un 87,7% de concordancia con el experto y un 91,9% de precisión sobre 2 registros de siesta diurna completos equivalentes a 6,3 horas, considerando características frecuenciales, morfológicas y criterios difusos. En este trabajo se presenta un detector que actúa sobre todo el registro de sueño, sin necesidad de establecer previamente el hipnograma, en que se aplican Descomposición Modal Empírica (Empirical Mode Descomposition, EMD),

1

JORNADAS CHILENAS DE INGENIERIA BIOMEDICA 2007 Transformada de Hilbert-Huang (THH), FFT y criterios expertos difusos. I.A. DESCOMPOSICIÓN MODAL EMPÍRICA (EMD)

LA EMD [11] permite descomponer de forma adaptiva señales no lineales y no estacionarias, en una colección de señales componentes denominadas Funciones Modales Intrínsecas (Intrinsic Mode Functions, IMFs). La EMD opera en forma iterativa, separando una serie temporal en una parte de alta frecuencia (IMF) y una parte con frecuencias menores o residuo. Luego el residuo se considera como la nueva serie temporal entrante, sobre la que se procede de la forma anterior. La serie temporal original puede ser reestablecida sumando los componentes que arroja la EMD.

En la figura 2 se observa un ejemplo de descomposición EMD sobre una señal EEG de 20 s en etapa SQ2. I.B. TRANSFORMADA DE HILBERT-HUANG

Una señal se representa usualmente en el dominio del tiempo o en el dominio de la frecuencia; la FFT es un ejemplo del segundo caso. En cambio, la Transforma de Hilbert (HT, por sus siglas en inglés) representa una señal en ambos dominios de forma simultánea. La HT de una señal x(t) se define como: HT ( x (t )) =

1

π

∞

x(τ ) dτ , t ( −τ ) −∞

∫

Señal EEG original en etapa SQ2, b)- i) Descomposición EMD en los modos IMF 1 a 8 de la señal EEG original, j) Residuo de la señal EEG.

El algoritmo EMD opera de la siguiente manera sobre una señal o serie de tiempo cualquiera: 1. Identificar la ventana de análisis. 2. Calcular los mínimos y máximos de la señal original (X(t)). 3. Interpolar sobre los mínimos obteniendo una envolvente inferior para la señal. 4. Interpolar sobre los máximos obteniendo una envolvente superior para la señal. 5. Calcular el promedio instantáneo entre las interpolaciones de mínimos y máximos, con lo que se construye la señal promedio (m1). 6. Calcular la diferencia instantánea entre la señal original y m1, con lo que se construye la primera señal residuo (r1 = X(t) - m1). 7. Tomar r1 como la nueva señal a procesar, y repetir los pasos 2 – 6, iterando hasta obtener una señal promedio de valor cero o al completar el límite de iteraciones (proceso de sifting). El residuo final resultante corresponde al primer modo IMF (IMF1). 8. Calcular la diferencia instantánea entre la señal original y el modo IMF obtenido. Tomar la nueva señal como la entrada para la construcción del siguiente modo IMF. 9. Repetir el proceso, iterando hasta obtener una señal residuo de valor cero o al completar el límite de modos.

(1)

que representa la convolución entre la señal x(t) y una función h(t)=1/πt. Esta operación introduce un desfase en el dominio frecuencial de -90º para frecuencias positivas y un desfase de 90º para frecuencias negativas. A partir de la HT se puede definir la señal analítica xsa(t) de x(t) correspondiente a la representación de una señal compleja a partir de una señal de valores reales: x sa (t ) = x(t ) + iHT ( x(t )). (2) Una señal analítica contiene sólo el espectro de Fourier positivo de la señal original. La señal analítica xsa(t) puede reescribirse como:

x sa (t ) = a(t )e iθ (t ) , con amplitud: a (t ) = x sa (t )

Fig. 2. Descomposición EMD para un trozo de registro EEG de 20 s. a)

(THH)

(3)

(4) y ángulo: θ (t) = arctan( HT ( x(t )) / x(t )), (5) donde a(t) corresponde a la amplitud instantánea de la señal x(t), y θ(t) representa la fase. Se puede definir la frecuencia instantánea w(t) como: dθ (t ) w(t ) = . (6) dt La THH [12] corresponde al proceso de calcular la HT de cada uno de los modos IMF y el residuo de la descomposición EMD para una señal x(t). La señal original x(t) puede ser reconstruida mediante N

i w (t ) dt x (t ) = ∑ a k (t ) e ∫ k , k =1

(7)

donde ak(t) corresponde a la frecuencia instantánea para el modo IMF k, y ∫ wk (t ) dt representa la fase del modo IMF k. II.

METODOLOGÍA

Los registros polisomnográficos fueron obtenidos utilizando un polígrafo Easy EEG-II (Cadwell, WA) de 32 canales, con el software versión 1.6.1. Los electrodos del EEG se posicionaron según el sistema internacional 10-20 para niños. Todos los canales fueron muestreados a 200 Hz. Para este trabajo se utilizaron las 2 derivaciones anteriores del EEG, F4C4, canal 1 (C1) y F3-C3, canal 3 (C3). Se utilizaron como conjunto de entrenamiento del sistema trozos de registros de distintos sujetos sanos de 10 años de edad en todos los estados y etapas de sueño, equivalentes

2

JORNADAS CHILENAS DE INGENIERIA BIOMEDICA 2007 aproximadamente a 39,7 horas. El conjunto de prueba usado estaba formado por 4 registros de sujetos distintos a los utilizados en el conjunto de entrenamiento, equivalentes aproximadamente a 34,2 horas de sueño nocturno. La validación de los resultados se realizó en base a la inspección visual de un experto médico del Laboratorio de Sueño del Instituto de Nutrición y Tecnología de los Alimentos (INTA) de la Universidad de Chile. El sistema de detección de HS se detalla en la figura 3.

Pbanda [a, b] i =

k ≤b

∑ P[k ]i ,

k ≥a

(9)

donde Pbanda[a,b] i es la potencia total en la banda [a,b] para la ventana de Hamming i, en que a y b son las frecuencias mínima y máxima de la banda respectivamente. 4. Avanzar la ventana de Hamming N/2 muestras y repetir pasos 1-3. 5. Detener proceso de cálculo de Pbanda[a,b] i al completar una zona de análisis EEG equivalente a 20 s de registro continuo. 6. Calcular la potencia promedio en la banda de interés en la zona de análisis EEG de 20 s como el promedio de las potencias totales de todas las ventanas contenidas en la zona de análisis de 20 s, según muestra la ecuación (10). − 1 nv P banda [a, b] = ∑ Pbanda [a, b]i , (10) nv i =1 donde nv corresponde al número de ventanas de Hamming contenidas en los 20 s de registro continuo. 7. Asociar a la zona de análisis EEG de 20 s el valor de −

Fig. 3. Diagrama de bloques del sistema de detección automática de patrones HS. Las derivaciones anteriores del EEG ingresan al Módulo 1 donde son procesadas utilizando la FFT para la detección de zonas de análisis en las cuales centrar la detección de HS. El Módulo 2 utiliza EMD-THH y criterios difusos y expertos sobre las zonas EEG de análisis definidas por el Módulo 1para generar candidatos a HS (HSc). El Módulo 3 filtra y valida la salida anterior basado en características extraídas de los HSc, entregando la detección final del sistema.

P banda [a, b] . 8. Repetir pasos 1-7 considerando un trozo de registro EEG de 20 s inmediatamente posterior. 9. Detener proceso hasta completar el registro. En una zona de análisis EEG de 20 s la potencia en la banda sigma debe superar un umbral y la potencia en la banda de alta frecuencia debe estar bajo otro umbral, las zonas de análisis que cumplen con los criterios anteriores ingresan al Módulo 2 para la detección de los HSc. La figura 4 muestra un ejemplo de detección de zonas de análisis EEG.

II.A. Módulo 1: Detección de zonas de análisis La definición de las zonas de análisis se realizó en base a un estudio del entorno del registro utilizando FFT sobre C1 y C3 del EEG. Se calculó la potencia espectral de las bandas sigma ([10, 15] Hz) y Alta Frecuencia ([30, 60] Hz), utilizando una ventana Hamming móvil de 2,56 s (512 muestras), con un traslape de 1,28 s (256 muestras) entre ventanas. El algoritmo de detección de zonas de análisis opera de la siguiente manera: 1. Calcular la FFT de la señal EEG en la ventana de Hamming i, obteniendo los coeficientes de la transformación. 2. Calcular la potencia espectral P[k]i como indica la ecuación (8). X [k ] i × X [k ]i * P[k ] i = , (8) N donde X[k] i corresponde a los coeficientes de la FFT para la ventana de Hamming i y N es el número de puntos de la ventana de Hamming i. 3. Calcular la potencia total en la banda de interés (sigma y alta frecuencia) como la sumatoria de los coeficientes de la potencia espectral (P[k]) contenidos en las bandas de interés, según la fórmula de la ecuación (9).

Fig. 4. El sistema de detección de zonas de análisis permite detectar las zonas del registro en las cuales debiera realizarse el análisis para la detección de HSc. a) Muestra el hipnograma para un registro, en el que se puede ver la evolución de los estados y etapas de sueño. b) Muestra la detección de de zonas para el C1del EEG, ‘1’ significa que la zona de análisis ingresa al Módulo 2, ‘0’ eliminar del análisis. Se puede observar que el sistema es capaz de detectar todas las etapas en las cuales aparecen HS (SQ2-4) y de eliminar del análisis gran parte de los estados y etapas en que no aparecen HS (vigilia, SQ1 y MOR).

II.B. Módulo 2: Generación candidatos Husos Sigma (HSc) Los HSc se generaron utilizando el algoritmo EMD-THH para obtener la amplitud instantánea (a(t)) y la frecuencia instantánea (w(t)) en las zonas de análisis definidas por el Módulo 1. Los cálculos de a(t) y w(t) se realizaron sobre los modos IMFs 2 y 3 de la EMD, ya que en estos modos se rescatan las componentes de frecuencia en el rango sigma. a(t) y w(t) son señales altamente ruidosas por lo que se debió 3

JORNADAS CHILENAS DE INGENIERIA BIOMEDICA 2007 aplicar un filtro pasabajos promediador de tipo rectangular de 0,125 s. A partir de a(t) y w(t) se crearon dos conceptos difusos, amplitud instantánea difusa (af(a)) y frecuencia instantánea difusa (wf(w)) utilizando los conjuntos difusos trapezoidales mostrados en la tabla I. TABLA I CONJUNTOS DIFUSOS PARA AMPLITUD Y FRECUENCIA INSTANTÁNEA Conjunto Difuso Variable Amplitud Instantánea Difusa (af(a)) Amplitud Instantánea

[10; 20; 200; 200]

Variable

Conjunto Difuso Frecuencia Instantánea Difusa (wf(w))

Frecuencia Instantánea

[9,5; 10,0; 15,0; 15,5]

con a en µV y w en Hz. Ambos conceptos difusos son aplicados sobre a(t) y w(t) de los modos IMFs 2 y 3. Se utilizó como norma-T el producto instantáneo (PI(t)) entre af(t) y wf(t). El proceso de generación de los HSc se realizó de la siguiente manera: 1. Las muestras de ambos modos, cuyos valores PI(t) sean menores a 0,25 son eliminadas. 2. De las muestras no eliminadas, aquellas cuya separación temporal sea menor a 0,05 s son unidas formando un pulso. Si la distancia entre muestras es mayor a 0,05 s la muestra es aislada y se comienza la generación de un nuevo pulso. 3. Los pulsos generados anteriormente son unidos si su distancia temporal es menor a 0,2 s formando un HSc. Si la distancia temporal entre pulsos es mayor a 0,2 s el pulso queda aislado. 4. Los HSc de duración menor a 0,5 s son eliminados. La figura 5 muestra un ejemplo del proceso de generación de los HSc sobre 20 s de registro EEG en etapa de sueño SQ2.

II.C. Módulo 3: Características Husos Sigma Los HSc generados son procesados para extraer características de tipo morfológicas y frecuenciales, para luego utilizar criterios de eliminación a partir de estas variables. El proceso de extracción se realiza a través de la detección de máximos y mínimos de la señal EEG original utilizando el método de mínimos cuadrados y cambios de pendiente. Las características extraídas a cada candidato son las siguientes: • Amplitud promedio HSc • Frecuencia promedio HSc • Duración HSc Este proceso permite el filtrado y validación de los HSc generados por el Módulo 2, de manera de refinar la detección. III.

Los resultados obtenidos en 4 registros de sueño nocturno completos (conjunto de prueba) se muestran a continuación. Se obtuvo un 90,6% de concordancia con la inspección visual del experto y un 74,1% de precisión en el conjunto de prueba. Los resultados muestran además, una baja dispersión tanto para la detección como la precisión del sistema. Los mejores resultados obtenidos muestran un 96,2% de concordancia con el experto y un 78,7% de precisión. Por otra parte los resultados más bajos entregados por el sistema muestran un 87,2% de concordancia en la detección y un 71,4% de precisión. Los resultados detallados se muestran en la tabla II. TABLA II RESULTADOS DE DETECCION DE HUSOS SIGMA PARA EL CONJUNTO DE PRUEBA Paciente Canal 1 2 3 4 TOTAL

1 3 1 3 1 3 1 3

Husos Sigma Concordancia ExpertoMarcas Detección Sistema Experto Sistema 2289 2617 2004 2059 2506 1796 1813 2322 1698 2161 2865 2079 1884 2155 1695 1860 2154 1642 1103 1394 998 1225 1579 1128 14394 17592 13040

IV.

Fig. 5. Proceso de generación de candidatos a Husos Sigma (HSc). a) muestra la señal original en etapa SQ2, donde se pueden observar 3 HS, b) Modo IMF 2 de la EMD de la señal original, c) w(t)para el modo IMF 2, d) a(t)para el modo IMF 2, e) Modo IMF 3 de la EMD de la señal original, f) w(t)para el modo IMF 3, g) a(t)para el modo IMF 3, h) corresponde a los HSc generados por el Módulo 2.

RESULTADOS

Marcadas, pero no detectadas 285 263 115 82 189 218 105 97 1354

Detectadas, pero no marcadas 613 710 624 786 460 512 396 451 4552

Porcentaje Detección

Precisión Sistema

87,5% 87,2% 93,7% 96,2% 90,0% 88,3% 90,5% 92,1% 90,6%

76,6% 71,7% 73,1% 72,6% 78,7% 76,2% 71,6% 71,4% 74,1%

DISCUSIÓN Y CONCLUSIONES

El sistema de detección de HS desarrollado muestra una concordancia con la inspección visual del experto (tasa de detección) de 90,6%, y una precisión de 74,1% en promedio, sobre el conjunto de prueba. Los resultados no son directamente comparables con nuestro trabajo previo [8], ya que los registros utilizados anteriormente corresponden a periodos de siesta diurna y no a sueño nocturno, por lo que la cantidad de HS son muy diferentes. No obstante lo anterior, cabe destacar que aunque los resultados no son sustancialmente mejores en comparación al trabajo previo, presentando un 2,9% de mejora en la detección, el nuevo sistema presenta ventajas principalmente en la localización de los HS a lo largo del registro, evitando de esta manera sobredimensionar las marcas generadas (principal problema del trabajo anterior). Otra diferencia importante entre el trabajo actual y el desarrollado en [8], radica en los conjuntos 4

JORNADAS CHILENAS DE INGENIERIA BIOMEDICA 2007 utilizados tanto para el entrenamiento como la prueba de los sistemas. En este trabajo el conjunto de entrenamiento equivale a 39,7 horas de registro de sueño nocturno, mientras que el conjunto de prueba equivale a cerca de 34,2 horas, en contraposición, el trabajo anterior fue entrenado utilizando un conjunto equivalente a cerca de 8 horas de registro y el conjunto de prueba equivalía a 6,3 horas. En relación al trabajo mostrado en [5], nuestro sistema presenta una mejoría en las tasas de detección, aumentando en un 20,6% la tasa de correctas detecciones, pero con una tasa de precisión menor, presentando una baja de 24,5%. La principal ventaja del sistema radica en la posibilidad de detectar HS sin requerir del hipnograma, a diferencia de otros trabajos ([5],[7]). Lo anterior se traduce en una disminución sustancial de los tiempos de experto requeridos, ya que para la generación de los hipnogramas se necesita un arduo trabajo. El sistema de detección será perfeccionado incorporando un algoritmo que permita seleccionar de manera automática los modos IMFs en el rango de frecuencias de interés. Por otra parte, se contempla la aplicación de clasificadores tipo redes neuronales y SVM.

for automated classification of sleep-waking states in infants”, Medical & Biological Engineering & Computing, vol. 40, pp. 105- 113, 2002. [11] G. Rilling, P. Flandrin, y P. Gonçalves, “On empirical mode decomposition and its algorithms”, en Proc. of the IEEE-EURASIP Workshop on nonlinear signal and image processing NSIP-03, vol. 1, 2003. [12] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, y H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis”, Proc. R. Soc. Lond. A, vol. 454, pp. 903–995, 1998.

AGRADECIMIENTOS Este trabajo fue financiado por CONICYT Chile, proyecto FONDECYT Nº 1070668, y proyecto NIH HD 33487. REFERENCIAS [1]

W. B. Martín, L. C. Johnson, S. S. Viglione, P. Naitoh, R. D. Joseph, J. D. Moses, “Pattern Recognition of EEG-EOG as a Technique for Allnight Sleep Stage Scoring”, Electroencephalogr. Clin. Neurophysiol., Vol. 32, pp. 417-427, 1972. [2] K. D. Nielsen, A. Kjaer, W. Jensen, T. Dyrby, L. Andreasen, J. Andersen, “Causal Probabilistic Network and Power Spectral Estimation Used in Sleep Stage Classification”, Methods Inf. Med., Vol. 36, pp. 345-348, 1997. [3] A. Rechtschaffen, A. Kales, “A Manual of Standardized Terminology, Techniques and Scoring System for Sleep Stages of Human Subjects”, UCLA, Brain Research Institute/ Brain Information Service, Los Angeles, USA, 1968. [4] S. Gais, M. Mölle, K. Helms, J. Born, “Learning-Dependent Increases in Sleep Spindle Density”, The Journal of Neuroscience, vol. 22, no. 15, pp. 6830-6834, 2002. [5] E. Huupponen, G. Gómez-Herrero, A. Saastamoinen, A. Värri, J. Hasan, S. Himanen, “Development and comparison of four sleep spindle detection methods”, Artificial Intelligence in Medicine, vol. 40, Issue 3, pp. 157-170. 2007. [6] E. Huupponen, A. Värri, S.-L. Himanen, J. Hasan, M. Lehtokangas, J. Saarinen, “Optimization of sigma amplitude threshold in sigma spindle detection”, J. Sleep Research, vol. 9, pp. 327-334, 2000. [7] D. Gorur, U. Halici, H. Aydin, G. Ongun, F. Ozgen, K. Leblebicioglu, “Sleep-spindles detection using short time Fourier transform and neural networks”, in Proc. IEEE International Joint Conference on Neural Networks, 2002. [8] Held C. M., L. Causa, P. Estévez, C. Pérez, M. Garrido, C. Algarín, P. Peirano, “Dual Approach for Automated Sleep Spindles Detection within EEG Background Activity in Infant Polysomnograms”, Proceedings of the 26th Annual International Conference of the IEEE EMBS, San Francisco, CA, USA, September 1-5, 2004. [9] C.M.Held, J.Heiss, P.A.Estevez, C.A.Perez, M.Garrido, C. Algarín, P.Peirano, "Extracting Fuzzy Rules From Polysomnographic Recordings For Infant Sleep Classification", IEEE Transactions on Biomedical Engineering, 0018-9294 Vol. 53 (10), pp. 1954 -1962, Oct 2006. [10] P. Estévez, C. M. Held, C. A. Holzmann, C. A. Perez, J. P. Pérez, J. Heiss, M. Garrido, P. Peirano, “Polysomnographic pattern recognition

5