CC BY
184
Изв. вузов «ПНД», т. 19, № 4, 2011
УДК 530.182(075.8)
ВЕИВЛЕТНЫИ АНАЛИЗ СОННЫХ ВЕРЕТЕН НА ЭЭГ И РАЗРАБОТКА МЕТОДА ИХ АВТОМАТИЧЕСКОЙ ДИАГНОСТИКИ*
В.В. Грубое, А.А. Овчинников, Е.Ю. Ситникова,
А.А. Короновский, А.Е. Храмов
Проведен детальный вейвлетный анализ электрической активности головного мозга, зарегистрированной во время сна у крыс с генетической предрасположенностью к абсанс-эпилепсии. Были использованы долговременные записи электроэнцефалограмм, где выявлены характерные особенности частотно-временной структуры сонных веретен - осцилляторных паттернов, являющихся электроэнцефалографическим коррелятом медленно-волнового сна. Работа проведена с использованием непрерывного вейвлетного преобразования. Предварительно были протестированы несколько наиболее распространенных вейвлетных базисов, среди которых комплексный вейвлет Морле оказался наиболее эффективным для выявления частотно-временных особенностей сонных веретен на электроэнцефалограммах. В дальнейшем вейвлет Морле был использован для создания системы автоматической диагностики сонных веретен на электроэнцефалограммах.
В результате удалось выделить два типа сонных веретен, имеющих сходную временную динамику, но различающихся по частотному составу. Выявлена сложная динамика основной частоты колебаний в течение сонного веретена. На основании полученных данных предложен метод автоматической диагностики сонных веретен, основанный на расчете величин мгновенных энергий вейвлетного преобразования в двух частотных диапазонах, характерных для двух типов сонных веретен. Тестирование метода показало высокую точность автоматической диагностики исследуемых событий на электроэнцефалограммах. Данный подход может быть использован в рутинных исследованиях электроэнцефалограмм, связанных с выделением и классификацией различных осцилляторных паттернов.
Ключевые слова: Непрерывное вейвлетное преобразование, ЭЭГ, эпилепсия, осциллятор-ный паттерн, сонные веретена, диагностика.
Введение
Современные математические методы анализа сложных колебательных процессов все чаще находят применение в самых различных областях естествознания, в том числе и в области физиологии и медицины [1-4]. Эти методы оказываются
‘Статья написана по материалам доклада на IX Международной школе «Хаотические автоколебания и образование структур», Саратов, Россия, 4-9 октября, 2010.
91
наиболее эффективными в исследованиях динамики активности нейронных сетей. Последние представляют собой крайне сложные объекты, состоящие из большого числа нейронов, где каждый нейрон можно рассматривать как отдельный колебательный элемент, характеризующийся собственной сложной динамикой [4-6]. Для получения информации о работе головного мозга традиционно используют методы регистрации электрической активности путем записи электроэнцефалограмм (ЭЭГ). ЭЭГ представляет собой усредненную сумму электрических токов, генерируемых группой нейронов, находящихся в области регистрирующего электрода [7, 8]. У человека регистрирующие электроды (сенсоры) обычно размещаются на поверхности головы. Для исследования животных существует практика вживления электродов в поверхностные и глубокие структуры головного мозга. Это дает возможность проведения более точных измерений, поскольку регистрирующие электроды находятся в непосредственном контакте с нервной тканью. Кроме того, запись ЭЭГ у животных можно проводить непрерывно на протяжении долгого периода времени (суток и более), что позволяет получить более полную информацию о мозговой активности.
Одна из задач электроэнцефалографических исследований связана с выявлением определенных форм ритмической активности. Традиционно на ЭЭГ выделяют несколько частотных диапазонов (альфа, бета, гамма, дельта и т.д.). Существуют четкие корреляции между появлением той или иной частотной составляющей на ЭЭГ и функциональным состоянием организма. Одним из наиболее многочисленных ос-цилляторных паттернов на ЭЭГ являются сонные веретена. Они наблюдаются во время второй фазы сна («медленный сон») и представляют собой колебания характерной «веретенообразной» формы (с увеличивающейся и далее уменьшающейся амплитудой осцилляций) с частотой 10...15 Гц и продолжительностью 0.5...1.5 с. Известно, что сонные веретена формируются вследствие синхронной работы нейронной сети, объединяющей нейроны коры и таламуса. Интерес к изучению сонных веретен продиктован их возможной связью с эпилептической активностью; в частности, существует гипотеза о трансформации сонных веретен в эпилептические пик-волновые разряды [9].
В настоящей работе была исследована частотно-временная структура сонных веретен на ЭЭГ у крыс с генетической формой абсанс-эпилепсии. По сути, на ЭЭГ у этих животных содержались осцилляторные паттерны эпилептической и неэпилептической природы. К первым относились пик-волновые разряды, характерные для абсанс-эпилепсии, ко вторым - сонные веретена.
Диагностика сонных веретен на ЭЭГ сопряжена с рядом технических сложностей из-за отсутствия строгих критериев их оценки на ЭЭГ, а также из-за существенного многообразия форм сонных веретен. Колебания электрического потенциала внутри отдельных веретен могут быть несимметричными (с преобладанием негативного компонента), могут появляться резкие всплески амплитуды (пики), может наблюдаться «дрейф» частоты основной ритмической составляющей веретена. Всё это порождает многообразие форм сонных веретен, создает дополнительные сложности на пути их анализа и затрудняет сопоставление данных, полученных в разных лабораториях. Наша работа направлена на разработку эффективных методов автоматической диагностики сонных веретен на ЭЭГ, а также принципов их стандартизации с использованием методов, развитых в нелинейной динамике. Представленные разработки могут найти применение в клинической электроэнцефалографии, в со-мнологии и фундаментальных исследованиях сна.
92
Основным инструментом в настоящей работе служило непрерывное вейвлет-ное преобразование [10-12]. Этот метод имеет ряд преимуществ и в настоящее время широко используется для анализа биомедицинских данных, например, в исследованиях реакции нервной системы на воздействие внешнего стимула [13, 14], ритмической активности респираторной и сердечно-сосудистой систем [15-18], синхронизации между различными областями головного мозга у больных эпилепсией [19, 20], для диагностики перемежающегося поведения в нейронных ансамблях [21-23], для выделения особенностей спайковой активности отдельных нейронов [24, 25]. На предварительном этапе настоящей работы было протестировано несколько базисных функций для непрерывного вейвлетного преобразования и выбран базис, дающий наибольшее частотно-временное разрешение применительно к сигналу ЭЭГ. Затем метод был оптимизирован для автоматического распознавания сонных веретен. В целом алгоритм распознавания веретен является модификацией предложенного нами ранее алгоритма выделения эпилептических пик-волновых разрядов [26-28].
1. Аппарат вейвлетного преобразования и исследуемые экспериментальные данные
Для выделения и анализа осцилляторных паттернов на ЭЭГ было использовано непрерывное вейвлетное преобразование, представляющее собой свертку исследуемого экспериментального сигнала x(t) и некоторой базисной функции [29]
W (s, т)
x(t)W*sT (t)dt
(1)
где * означает комплексное сопряжение.
Базисная функция tys,T(t) может быть получена из материнского вейвлета с помощью следующего преобразования:
^S,T (t)
7 ф0
t — т s
(2)
где s - временной масштаб, определяющий растяжение или сжатие материнской функции, т - временной сдвиг вейвлетного преобразования и фо - материнский вейвлет. Часто вместо временного масштаба s рассматривают частоту fs = 1/s, которую вводят по аналогии с частотой спектра Фурье. В то же самое время, следует обратить внимание на то, что интуитивно ожидаемое соответствие fs = f между величиной fs вейвлетного преобразования и частотой f преобразования Фурье, в общем случае не справедливо и при трактовании полученных результатов это следует учитывать [29].
Для выявления особенностей частотной динамики сигнала в данный момент времени tфикс удобно использовать мгновенные распределения амплитуды сигнала по вейвлетному спектру
Efs = \W (fs,t = tфикс) |. (3)
В работе были использованы данные ЭЭГ, зарегистрированные в коре и в таламусе у восьми крыс эпилептической линии WAG/Rij [30]. Для анализа были
93
Рис. 1. Пример типичной записи многоканальной эпилептической сонной ЭЭГ, использованной в настоящем исследовании, а также схема расположения регистрирующих электродов
отобраны фрагменты ЭЭГ во время сна с хорошо выраженными сонными веретенами. Характерные записи многоканальной «сонной» ЭЭГ приведены на рис. 1, а. Расположение регистрирующих электродов представлено схематически на рис. 1, б. Исследуемые паттерны - сонные веретена - отмечены на временных реализациях штриховыми рамками. Хорошо видно, что наиболее ярко сонные веретена выделяются в лобной коре (канал Fr), поэтому далее все исследования проводились со скалярным сигналом, снимаемым именно из этой области головного мозга.
Первоначально были протестированы несколько материнских базисных функций для непрерывного вейвлетного преобразования «сонной» ЭЭГ и выявлены достоинства и недостатки каждой из них. Мы сосредоточили внимание на анализе ЭЭГ с помощью действительного - MHAT [31] - и комплексных - Паула [32] и Морле [33] - материнских вейвлетов, поскольку эти вейвлетные функции наиболее часто используются при исследованиях структуры сигналов. Результаты этого анализа представлены в разделе 2. Необходимо иметь в виду, что каждый вейвлетный базис характеризуется собственным соотношением между частотами fs вейвлетного преобразования и фурье-частотами, с которыми исследователь имеет дело на практике. Соответственно для MHAT-вейвлета данное соотношение имеет вид fs ~ 3.97f, для вейвлета Паула с коэффициентом m = 4 - fs w 0.71f, и для вейвлета Морле с собственной частотой Q = 2п - fs f [29].
На исследованных записях ЭЭГ была проведена «разметка», то есть выделены моменты начала и окончания сонных веретен и пик-волновых разрядов (эпилептических событий). В дальнейшем была проведена оценка эффективности автоматического распознавания сонных веретен путем сравнения экспертной разметки с результатами вейвлетного анализа.
2. Вейвлетный анализ «сонной» ЭЭГ
Результаты вейвлетного анализа ЭЭГ с большим количеством сонных веретен, зарегистрированной в лобной коре во время сна, показаны на рис. 2. Здесь представлен характерный отрезок ЭЭГ и соответствующие амплитудные вейвлет-ные спектры, полученные с помощью материнских MHAT (а), Паула (б) и Морле (в) вейвлетных базисов. Штриховые рамки обозначают сонные веретена на ЭЭГ и соответствующие им участки вейвлетных спектров. Учитывая, что частоты сонных
94
г о 4 8 12 t,c
Рис. 2. Отрезок сонной ЭЭГ из фронтальной области (а) и вейвлетные спектры, построенные с различными материнскими вейвлетами: MHAT (б), Паула (в) и Морле (г)
95
веретен лежат в диапазоне 10...15 Гц (альфа частоты), этот частотной диапазон был выбран для анализа распределения коэффициентов вейвлетного преобразования.
Анализ и сравнение различных типов материнских базисов показал, что оптимальным для выявления особенностей частотно-временной структуры сонных веретен на ЭЭГ является вейвлет Морле, на основании чего было сделано заключение об использовании именно этого базиса для создания системы автоматического выделения сонных веретен. Действительно, если рассмотреть вейвлетный спектр, полученный с помощью MHAT-вейвлета, мы обнаружим, что частотное разрешение данного вейвлета невелико, и диапазон в частотной области сонных веретен оказывается достаточно широким и перекрывается с диапазонами, соответствующими другим осцилляторным событиям на ЭЭГ. В то же время, особенности действительного MHAT-вейвлета таковы, что резкие пики и прочие события на ЭЭГ демонстрируют на вейвлетном спектре мощные всплески амплитуды вейвлетных коэффициентов в интересующем нас диапазоне частот. Кроме того, форма вейвлетного спектра сонных веретен не позволяет отследить особенности динамики их частоты и амплитуды. К достоинствам MHAT-вейвлета стоит отнести хорошее временное разрешение, которое определяет резкое увеличение амплитуды коэффициентов вейвлетного спектра при возникновении соответствующего ритма на ЭЭГ. Однако из-за плохого частотного разрешения сложно классифицировать резкие скачки амплитуды вейвлетных спектров.
Аналогичная ситуация наблюдается и в случае с комплексным материнским вейвлетом Паула. Его частотное разрешение также не дает возможности точно определить особенности динамики сигнала в области альфа частот. В рассматриваемом частотном диапазоне на вейвлетной поверхности наблюдается много высокоамплитудных составляющих, не являющихся сонными веретенами, что затрудняет проведение анализа.
Иную картину демонстрирует вейвлет Морле, который сочетает в себе хорошее разрешение как в частотной, так и во временной областях сигнала. В области вейвлетного спектра, соответствующей альфа-диапазону, можно легко локализовать сонные веретена как на шкале времени, так и на шкале частот. Это дает возможность отследить на ЭЭГ изменения частоты и амплитуды колебаний, характерных для сонного веретена, в автоматическом режиме.
Проведенный анализ большого числа вейвлетных спектров сонных веретен показал, что на записях ЭЭГ у крыс эпилептической линии выделяется несколько характерных типов сонных веретен по критерию их частотной динамики, а именно основной частоты вейвлетного спектра. Как уже отмечалось во введении, сонные веретена имеют частоту порядка 10...15 Гц, при этом их форма весьма разнообразна, поэтому базовые частоты существенно варьируются в данном частотном диапазоне от веретена к веретену. Однако на вейвлетных спектрах, построенных с помощью материнского вейвлета Морле, были обнаружены события, форма которых была близка к форме сонных веретен, но базовые частоты лежали в области, характерной для эпилептических пик-волновых разрядов (7...12 Гц, см., например [34]). Данный факт удобно проиллюстрировать мгновенными вейвлетными спектрами E (fs) типичных сонных веретен.
На рис. 3 показаны характерные распределения амплитуды вейвлетных коэффициентов, а также осциллограммы типичных паттернов ЭЭГ, которые использова-
96
Рис. 3. Распределения амплитуды коэффициентов вейвлетного преобразования по частотам для различных отрезков ЭЭГ, соответствующих фоновой активности (1), типичному (2) и атипичному (проэпилептическому) (3) сонным веретенам, построенные с различными материнскими вейвлетами: MHAT (а), Паула (б) и Морле (в)
лись для построения мгновенных вейвлетных спектров. Представлены результаты вейвлетных преобразований с базовыми вейвлетными функциями MHAT (а), Паула (б) и Морле (в). Первая колонка на рис. 3 соответствует фоновой ЭЭГ (состояние десинхронизации ЭЭГ). Хорошо видно, что в это время на ЭЭГ наблюдаются колебания в широком диапазоне частот, но с очень малой амплитудой. Сравнительно большой амплитудой вейветного спектра характеризуются только колебания с частотами, соответствующими диапазону 0...5 Гц. Вторая колонка на рис. 3 иллюстрирует спектр типичного сонного веретена. Видно, что максимум амплитуды вейвлетного спектра приходится на частоты 10...12 Гц, что соответствует классическому определению частотного состава сонного веретена. Следует отметить значительно большую (по сравнению с фоновой ЭЭГ) энергию, приходящуюся на спектр сонного веретена, что находит отражение в больших значениях мгновенного спектра E. В третьей колонке на рис. 3 представлен типичный осцилляторный паттерн, имеющий форму сонного веретена, но более низкую частоту (в диапазоне 6...10 Гц), соответствующую частоте эпилептического пик-волнового разряда. Подобные события на ЭЭГ
97
Таблица 1
№ живот-ного fs Гц Средняя длительность, с Wmax W) N УЛкр ^2кр
Fs! Fs2 Fsl Fs2 Fsl Fs2
1 8.38 ± 1.60 13.64 ± 1.70 0.63 ± 0.25 0.52 ± 0.44 1.35 ± 0.80 0.40 0.40 0.18 0.06 3.63
2 7.56 ± 1.18 14.3 ± 1.70 0.59 ± 0.27 0.22 ± 0.11 0.22 ± 0.12 0.20 0.20 0.10 0.06 1.99
3 8.17 ± 0.93 13.3 ± 1.26 0.47 ± 0.12 0.27 ± 0.15 0.25 ± 0.07 0.25 0.20 0.08 0.05 1.96
4 8.07 ± 0.73 15 ± 1.00 1.23 ± 0.73 0.87 ± 0.76 0.47 ± 0.21 0.40 0.40 0.15 0.10 1.44
были названы «атипичными» сонными веретенами. Максимум вейвлетного спектра «атипичного» веретена находился в области частот, близких к 8 Гц. Усредненные величины максимальных значений коэффициентов вейвлетного преобразования Wmax в каждом частотном диапазоне представлены в табл. 1.
Характерные особенности каждого из четырех анализируемых животных, результаты вейлетного анализа частотно-временных особенностей сонных веретен с использованием вейвлета Морле, а также параметры описанного в разделе 3 метода приведены в табл. 1. В таблице представлены следующие данные:
• доминирующие частоты fs колебаний сонного веретена в диапазонах частот Fs\ Е [5...9] Гц и Fs2 Е [10...16] Гц, которые, как будет показано ниже, имеют наибольшую мощность на вейвлетном спектре;
• средняя длительность сонного веретена каждого животного;
• усредненные величины максимальных значений Wmax вейвлетных коэффициентов в каждом характерном частотном диапазоне в течение сонного веретена;
• средние величины (W^,) вейвлетных коэффициентов в каждом характерном частотном диапазоне фоновой ЭЭГ;
• среднее число N сонных веретен, наблюдающихся на ЭЭГ во время сна на интервале 10 с;
• параметры метода автоматической диагностики сонных веретен (см. раздел 3 статьи).
Сравнивая мгновенные вейвлетные спектры, построенные с помощью различных вейвлетных базисов, можно сделать вывод, что только вейвлет Морле дает возможность выделить и различить два обнаруженных типа сонных веретен. Хорошо видно, что использование MHAT и комплексного вейвлета Паула в качестве базисных функций не обеспечивает необходимого разрешения по частоте (большая ширина соответствующего максимума в спектре на рис. 3, а и б), поэтому два пика на частотах 12 и 8 Гц практически сливаются и оказываются неразличимыми с точки зрения вейвлетного анализа. Следует также отметить, что энергия «атипичного» веретена оказывается существенно (в 2-3 раза) меньшей, чем энергия типичного сонного веретена; из-за этого «атипичные» веретена сложно выделить на фоне десинхронизованной («фоновой») ЭЭГ и отличить от более часто встречающихся типичных веретен. Данный факт также значительно усложняет автоматическое выделение сонных веретен с использованием вышеназванных типов материнских вейвлетов.
В то же самое время вейвлетный анализ с материнским комплексным вейвлетом Морле позволяет эффективно выделить два типа сонных веретен: во-первых,
98
благодаря хорошему частотному разрешению, пики, соответствующие двум типам веретен, на вейвлетных спектрах не перекрываются, во-вторых, нормировка амплитуды вейвлетного спектра такова, что амплитуды обоих типов веретен оказываются сравнимы друг с другом. Это делает возможным решение задачи о распознавании сонных веретен на ЭЭГ с использованием материнского вейвлета Морле.
Следует отметить, что сбалансированное частотно-временное разрешение, обеспечиваемое использованием комплексного вейвлета Морле с Q = 2п, позволяет провести эффективный анализ частотной динамики в течение сонного веретена, что было бы невозможно сделать с использованием других вейвлетных базисов. На рис. 4 показан отрезок ЭЭГ (а), соответствующий вейвлетный спектр (б), полученный с помощью комплексного вейвлета Морле, и мгновенные распределения амплитуды (в) вейвлетного спектра, соответствующие моментам времени 1-3, выделенным на рис. 4, а стрелками. На приведенном отрезке можно выделить два сонных
Рис. 4. Динамика частотного состава в течение сонного веретена: временная реализация из фронтальной области (а), вейвлетный спектр (б) и распределение амплитуды коэффициентов вейвлетного преобразования в различные моменты времени (отмечены стрелками на рис. а) (в)
99
веретена, но первое, появляющееся в момент времени t = 4 c, является коротким и слабо выраженным на фоновой ЭЭГ, поэтому анализ проводился для второго сонного веретена, возникающего в момент времени t = 9... 11 c и отмеченного стрелками. Наиболее важным результатом такого детального анализа является тот факт, что основная частота в течение сонного веретена сильно меняется, увеличиваясь к концу веретена. Это хорошо можно проследить на вейвлетной поверхности и на представленных мгновенных спектрах. Видно, что колебания в начале сонного веретена наблюдаются на частоте порядка 7.5 Гц (кривая 1 на рис. 4, в), затем растет амплитуда колебаний и наблюдается медленный сдвиг частоты в сторону больших частот (кривая 2). Одновременно с этим в спектре имеют место более быстрые колебания с меньшей амплитудой и частотой порядка 14 Гц. Однако к концу сонного веретена происходит быстрое изменение (увеличение) основной частоты до 12 Гц, в то же самое время, более высокочастотная динамика перестает наблюдаться, однако ритм 7.5 Гц, который был в самом начале сонного веретена, остается, хотя и с существенно меньшей амплитудой.
Рост основной частоты (в ряде случаев быстрый, как в рассматриваемом примере, в ряде случаев - более монотонный и плавный) оказывается типичным для сонных веретен на эпилептической ЭЭГ. Это отличает динамику сонных веретен от пик-волновых разрядов, которые также характеризуются сильными изменениями частоты внутри разряда. Однако частота пик-волнового разряда уменьшается, а не увеличивается к концу последнего, как в сонных веретенах. Следует также отметить, что такая сложная перестройка частоты колебаний в течение сонного веретена может существенно усложнить выделение типичного частотно-временного паттерна сонного веретена и, как следствие, создать определенные трудности при автоматической диагностике и выделении сонных веретен. Как было отмечено в работе [34], это связано с большей вариабельностью сонных веретен по сравнению с другими типичными осцилляторными паттернами на ЭЭГ, в частности, пик-волновыми разрядами [35].
3. Автоматическое выделение сонных веретен на ЭЭГ на основе анализа распределения энергии вейвлетного спектра
Рассмотрим возможность построения автоматической экспертной системы, которая позволила бы эффективно и с высокой точностью выделять сонные веретена на записях ЭЭГ на основе описанных выше частотно-временных особенностей сонных веретен. Известен ряд оригинальных методов анализа динамики нейронных ансамблей на базе ЭЭГ с помощью вейвлетного преобразования. Так, в работе [36] описывается метод классификации ЭЭГ с помощью оценки значений вейвлетпакетных коэффициентов в определенных частотных диапазонах. В работе [37] вей-влетное преобразование было применено для предварительной подготовки сигнала, а именно, для очистки от артефактов посредством дискриминации малозначимых вейвлет-коэффициентов, полученных на основе разложения сигнала по ортогональным вейвлет-пакетам. В работе [38] предложен метод анализа переходных процессов, наблюдаемых на ЭЭГ при фотостимуляции, основанный на наблюдении за динамикой спектральных компонент на вейвлетной поверхности. Широко применяет-
100
ся вейвлетный анализ и для выделения сигналов отдельных нейронов при записи электрической активности малых групп нейронов [24, 25]. Важно отметить, что наиболее эффективным способом в данном контексте является анализ распределения энергии непрерывного вейвлетного преобразования по частотам и выделение критериев, которые позволили бы описать характерные частотно-временные образы того или иного осцилляторного паттерна [26, 34]. Подобные методы оказываются эффективными для обработки уже записанных сигналов, а также позволяют детектировать осцилляторные паттерны в режиме реального времени (см., например, [27, 28]).
Алгоритм распознавания сонных веретен, использованный в рамках настоящего исследования, был впервые использован в работах [26, 34] для автоматической диагностики пик-волновых разрядов. Как уже отмечалось в разделе 2 статьи, вей-влетные спектры сонных веретен характеризуются пиком мощности в диапазоне частот 6...15 Гц, при этом базовая частота в спектре сонного веретена имеет тенденцию увеличиваться в пределах данного частотного диапазона от начала к концу сонного веретена. Основная идея использованных ранее методов [26, 34] заключалась в том, что исследуемый временной ряд (ЭЭГ) подвергался вейвлетному преобразованию, затем проводился расчет мгновенной энергии преобразования w(t) в характерном частотном диапазоне Fs:
w(t) = j\Ws,t\ds. (4)
Fs
В случае, когда значение w(t) превышало некоторый критический уровень Wkp (w(t) > w^p), делался вывод о наличии искомого паттерна на ЭЭГ.
В работах [26, 34] для выделения пик-волнового разряда был использован частотный диапазон Fs £ [35...100] Гц. В настоящей работе выделение сонных веретен проводили с учетом особенностей их частотного состава, вследствие чего были использованы два частотных диапазона Fsi £ [5...9] Гц и Fs2 £ [10...16] Гц, что позволило детектировать как обычные, так и «атипичные» сонные веретена (см. распределения энергии вейвлетного спектра, рис. 3, в). Вывод о наличии на ЭЭГ сонного веретена делался на основе анализа энергий w1(t) и w2(t) в каждом из диапазонов Fsi и Fs2, когда величины wi(t) и w2(t) превышали некоторые пороговые значения w^ и w^ (для первого и второго диапазонов, соответственно)
wi(t) > wi^ V W2(t) > W2^. (5)
Отличительной чертой вейвлетного преобразования сигнала ЭЭГ, содержащего сонное веретено, являлась сильная изрезанность функции энергии вейвлетного преобразования w(t) (4) от времени. Это обусловлено обнаруженными изменениями частоты в течение сонного веретена и вариациями амплитуды колебаний. Поэтому для определения момента прекращения сонного веретена был использован следующий критерий:
wi(t) < W0 Л W2(t) < wo, (6)
где wo - средний уровень мощности вейвлетного спектра, рассчитанный для фоновой ЭЭГ в исследуемом частотном диапазоне. Данный подход позволил избежать ложного обнаружения нескольких «коротких» сонных веретен, когда на самом деле существовало одно длинное веретено.
101
Наконец, необходимо отметить, что ЭЭГ является сложным сигналом, в котором могут появляться отдельные всплески активности в различных диапазонах частот, что отражает динамику функционирования нейронных ансамблей коры головного мозга (в частности, К-комплексы, альфа-подобная активность, комплексы дельта волн). Подобные события могут вызывать резкое кратковременное увеличение мгновенной энергии преобразования, в том числе и в рассматриваемых частотных диапазонах Fs\ и Fs2, и, соответственно, вызывать ложное детектирование сонного веретена. Поскольку подобные всплески энергии представляются достаточно частым событием, в вышеописанный алгоритм был внесен ряд усовершенствований для повышения точности анализа. Так, при диагностике колебаний в соотношениях (5) и (6) при сравнении величин w\ и W2 с пороговыми значениями wi^ и ^2кр использовали не мгновенные значения энергии преобразования (4), а значения, усредненные по некоторому интервалу времени:
Численный анализ показал, что оптимальной шириной окна, позволяющей снизить влияние коротких артефактов на качество распознавания веретен, является величина
Рис. 5 наглядно иллюстрирует принципы работы описанного метода (4)-(7) и результаты выделения сонных веретен. На рисунке представлен типичный отрезок «сонной» ЭЭГ с несколькими веретенами (а), соответствующий вейвлетный спектр (б) и энергия вейвлетного преобразования w1(t) и w2(t) в диапазонах Fs1 и Fs2 (в), а также сонные веретена, выделенные экспертом (г) и автоматически (д). Горизонтальными линиями на рис. 5, г показаны пороговые значения энергий wi^ и w2Kp, которые были определены эмпирически. Видно, что метод, основанный на анализе энергии ЭЭГ в двух частотных диапазонах, характерных для сонно-веретенной активности, обладает весьма высокой точностью, сопоставимой с экспертной разметкой. В указанном временном интервале наблюдалась единственная ошибка детектирования в момент времени t = 257 с, когда метод демонстрирует наличие сонного веретена, в то время как эксперт не отметил данное событие как веретено. Важным аспектом предложенного метода является то, что критические пороги энергий для двух рассматриваемых частотных диапазонов могут различаться (или быть одинаковыми), что позволяет гибко настроить алгоритм для максимальной точности выделения интересующих событий.
Оптимально подобранные параметры метода для каждого из животных представлены в табл. 1. Именно эти параметры были использованы далее для автоматизированной обработки записей ЭЭГ, содержащих сонные веретена.
Вычисления проводились на ЭВМ Pentium Quad 2.4 ГГц, 2 Гб ОЗУ, ОС Windows XP SP3. Программа анализа была разработана в свободно распространяемой среде программирования Python 2.7. Для обработки записи ЭЭГ с длительностью 2 часа требовалось около 4-х часов машинного времени. Так как расчет для продолжительной записи ЭЭГ занимает значительное машинное время, то представляется целесообразным разбиение одной длинной записи ЭЭГ на несколько более коротких отрезков, которые могут обрабатываться параллельно на различных ядрах процессора или на различных компьютерах с использованием, например, Grid-технологии. Это позволяет значительно снизить затраты времени, требуемого для обработки одной анализируемой записи.
(7)
т
T = 0.5 с.
102
Рис. 5. Отрезок сонной ЭЭГ из фронтальной области (а), вейвлетный спектр (б), мгновенные энергии вейвлетного преобразования в диапазонах Fs1 (штриховая линия) и Fs2 (сплошная линия) (в), выделенные экспертом (г) и предложенным в работе методом (д) сонные веретена. Горизонтальными линиями на (в) показаны пороги энергий ад1кр (штриховая линия) и ад2кр (сплошная линия)
На основании собранной статистики по качеству распознавания сонных веретен на ЭЭГ анализировались такие важные статистические характеристики автоматической диагностики сонных веретен на ЭЭГ, как уровень значимости в и мощность критерия 8 [39, 40]:
8
Ntp
(Ntp + Nfn )
х 100%,
в
Ntp
(Ntp + Nfp )
х 100%,
(8)
где Ntp - число верно распознанных событий, Nfp - число ложно распознанных событий, то есть число событий, которые определялись автоматически как сонное веретено, но не были распознаны экспертом, и Nfn - число пропущенных сонных веретен. Данные статистические параметры часто используются в качестве характеристик принятия «бинарного» решения (в данном случае наличие/отсутствие сонного веретена) на базе некоего критерия, который с некоторой вероятностью может давать ложный результат. Первая характеристика (8) позволяет оценить чувствительность метода, то есть процент распознанных сонных веретен по отношению к
103
Таблица 2
№ животного NTV Nfp Nfp 8,% в,%
1 38 6 1 97 86
2 39 12 30 76 55
3 30 3 12 90 64
4 27 2 4 93 87
Среднее 33.5 ± 5.9 6.3 ± 4.0 11.3 ± 4.0 89 ± 9.1 73 ± 16.0
общему числу сонных веретен на ЭЭГ; вторая (в) - процент событий, верно определенных как сонные веретена от общего числа событий, диагностированных как сонные веретена.
Результаты статистической проверки метода представлены в табл. 2, где видно, что уровень значимости 8 метода был порядка 90%, среднее значение критерия мощности в ~ 73%. Принимая во внимание сильную вариабельность и разнообразие паттерна сонного веретена на ЭЭГ, можно считать, что эффективность метода оказывается достаточно большой. Наихудшие результаты использования автоматического распознавания сонных веретен были получены у крысы № 2. ЭЭГ данной крысы характеризовалась очень низкой амплитудой, и сонные веретена были слабо различимы на фоновой записи даже экспертом. Можно предположить, что в данном случае пришлось иметь дело с примером неудачной записи ЭЭГ, что могло быть связано с индивидуальными особенностями животного, дефектом электродов или недостаточностью контакта c электродами.
Заключение
В представленной работе был проведен анализ электрической активности головного мозга, зарегистрированной во время сна у крыс с генетической предрасположенностью к абсанс-эпилепсии. Были исследованы характерные особенности частотно-временной структуры сонных веретен - особых осцилляторных паттернов, являющихся электроэнцефалографическим коррелятом медленно-волнового сна.
В ходе разработки технических средств для обработки сигнала ЭЭГ с помощью непрерывного вейвлетного преобразования было проведено тестирование различных материнских функций. Показано, что комплексный вейвлет Морле является наиболее эффективным базисным вейвлетом, позволившим выявить характерную частотно-временную структуру сонных веретен на ЭЭГ. Данный вывод хорошо согласуется с предыдущими результатами наших работ [26-28,34,41] по изучению других паттернов на ЭЭГ животных, больных абсанс-эпилепсией. Как следствие, вейвлет Морле был использован для создания системы автоматической диагностики сонных веретен на ЭЭГ.
Были обнаружены два типа сонных веретен, имевших сходную временную динамику, но разный частотный состав. Выявлена сложная динамика основной частоты в течение сонного веретена.
Основываясь на результатах работ [23,26,27,34], был разработан новый метод автоматической диагностики сонных веретен. Метод базируется на результатах непрерывного вейвлетного преобразования анализируемых сигналов и предусматри-
104
вает расчет величин мгновенных энергий вейвлетного преобразования в двух частотных диапазонах, характерных для выделяемых осцилляторных паттернов. Данный метод намного проще в реализации, чем метод, предложенный ранее в работе [34] для выделения сонных веретен и основанный на построении адаптивных вейвлетных базисов. Тем не менее, несмотря на простоту, предложенный нами метод продемонстрировал высокую точность и селективность в отношении сонных веретен.
Следует отметить, основываясь на результатах работы [27], посвященной выделению пик-волновых эпилептических разрядов на ЭЭГ в реальном времени с помощью алгоритма, сходного с рассмотренным в данной работе, что предложенный метод может быть в перспективе адаптирован и использован для диагностики сонных веретен на ЭЭГ во время записи электрической активности головного мозга в режиме реального времени.
Работа поддержана Президентской программой поддержки ведущих научных школ (проект НШ-3407.2010.2), а также ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы.
Библиографический список
1. Абарбанель Г.Д.И., Рабинович М.И., Селверстон А. и др. Синхронизация в нейронных ансамблях // Успехи физических наук. 1996. Т. 166. С. 363.
2. Mosekilde E., Maistrenko Yu., Postnov D.E. Chaotic synchronization, applications to living systems. Singapore: World Scientific. 2002.
3. Безручко Б.П., Пономаренко В.И., Прохоров и др. Моделирование и диагностика взаимодействия нелинейных колебательных систем по хаотическим временным рядам (приложения в нейрофизиологии) // Успехи физических наук. 2008. Т. 178. С. 323.
4. Некоркин В.И. Нелинейные колебания и волны в нейродинамике // Успехи физических наук. 2008. Т. 178. С. 313.
5. Steriade M., Deschenes M. The thalamus as a neuronal oscillator // Brain Res. Rev. 1984. Vol. 8. P. 1.
6. Destexhe A., Sejnowski T/.Thalamocortical assemblies: How ion channels, single neurons and large-scale networks organize sleep oscillations. Oxford University Press, 2001.
7. Niedermeyer E., Silva F.L. Electroencephalography: Basic principles, clinical applications, and related fields. Lippincot Williams & Wilkins. 2004.
8. Nunez P.L., Srinivasan R. Electric fields of the brain: The neurophysics of EEG. Oxford University Press, 1981.
9. Kostopoulos G.K. Spike-and-wave discharges of absence seizures as a transformation of sleep spindles: the continuing development of a hypothesis // Clinical Neurophysiology. 2000. Vol. 111. Suppl 2. S27-38.
10. Holschneider M. Wavelets: An analysis tool. Oxford University Press, 1995.
11. Aldroubi A., Unser M. Wavelets in medicine and biology. CRC-Press, 1996.
12. Daubechies I. Orthonormal bases of compactly supported wavelets // Communications on Pure and Applied Mathematics. 1988. Vol. 41. P. 909.
13. Tass P.A. et al. Detection of n:m phase locking from noisy data: Application to magnetoencephalography // Phys. Rev. Lett. 1998. Vol. 81. P. 3291.
105
14. Tass P.A., Fieseler T., Dammers J. et al. Synchronization tomography: A method for three-dimensional localization of phase synchronized neuronal populations in the human brain using magnetoencephalography // Phys. Rev. Lett. 2003. Vol. 90, P. 088101.
15. Anishchenko V.S., Balanov A.G., Janson N.B. et al. Entrainment between heart rate and weak noninvasive forcing // Int. J. Bifurcation and Chaos. 2000. Vol. 10. P. 2339.
16. Prokhorov M.D., Ponomarenko V.I., Gridnev V.I. et al. Synchronization between main rhythmic processes in the human cardiovascular system // Phys. Rev. E. 2003. Vol. 68. P. 041913.
17. Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. Detecting synchronization of self-sustained oscillators by external driving with varying frequency // Phys. Rev. E. 2006. Vol. 73. P. 026208.
18. Короновский А.А., Пономаренко В.И., Прохоров М.Д., Храмов А.Е. Метод исследования синхронизации автоколебаний по унивариантным данным с использованием непрерывного вейвлетного анализа // ЖТФ. 2007. T. 77, № 9. С. 6.
19. Meinecke F.C., Ziehe A., Kurths J., Muller K.-R. Measuring phase synchronization of superimposed signals // Phys. Rev. Lett. 2005. Vol. 94. P. 084102.
20. Chavez M., Adam C., Navarro V et al. On the intrinsic time scales involved in synchronization: A data-driven approach // Chaos. 2005. Vol. 15. P. 023904.
21. Velazquez J.L.P., Khosravani H., Lozano A. et al. Type III intermittency in human partial epilepsy // European Journal of Neuroscience. 1999. Vol. 11. P. 2571.
22. Короновский А.А., Кузнецова Г.Д., Мидзяновская И.С., Ситникова Е.Ю., Тру-бецков Д.И., Храмов А.Е. Закономерности перемежающегося поведения в спонтанной неконвульсивной судорожной активности у крыс // ДАН. 2006. T. 409.
C. 274.
23. Hramov A.E., Koronovskii A.A., Midzyanovskaya I.S. et al. On-off intermittency in time series of spontaneous paroxysmal activity in rats with genetic absence epilepsy // Chaos. 2006. Vol. 16. P. 043111.
24. Sosnovtseva O.V., Pavlov A.N., Mosekilde E., Yip K.-P., Holstein-Rathlou N.-H., Marsh D.J. Synchronization among mechanisms of renal autoregulation is reduced in hypertensive rats // American Journal of Physiology (Renal Physiology). 2007. Vol. 293. P. F1545.
25. Sosnovtseva O.V., Pavlov A.N., Mosekilde E., Holstein-Rathlou N.-H. Synchronization phenomena in multimode dynamics of coupled nephrons // Известия вузов. Прикладная нелинейная динамика. 2003. T. 11, №3. C. 133.
26. Короновский А.А., Минюхин И.М., Тыщенко А.А., Храмов А.Е., Мидзяновская И.С., Ситникова Е.Ю., ван Луйтелаар Ж., ван Рижн С.М. Применение непрерывного вейвлет-преобразования для анализа перемежающегося поведения // Изв. вузов. Прикладная нелинейная динамика. 2007. T. 15, №4. C. 34.
27. Ovchinnikov A.A., Luttjohann A., Hramov A.E., van Luijtelaar G. An algorithm for real-time detection of spike-wave discharges in rodents // Journal of Neuroscience Methods. 2010. Vol. 194. P. 172.
28. Овчинников А.А., Храмов А.Е., Люттьеханн А., Короновский А.А., ван Луйте-лаар Ж. Метод диагностики характерных паттернов на наблюдаемых временных рядах и его экспериментальная реализация в режиме реального времени применительно к нейрофизиологическим сигналам // ЖТФ. 2011. T. 81. C. 3.
29. Короновский А.А., Храмов А.Е. Непрерывный вейвлетный анализ и его приложения. М.: Физматлит, 2003.
106
30. Coenen A.M., van Luijtelaar E.L. Pharmacological dissociation of EEG and behavior: A basic problem in sleep-wake classification // Sleep. 1991. Vol. 14. P. 464.
31. Mallat S.G. Multiresolution approximations and wavelets orthonormal bases of L2(R) // Trans. Amer. Soc. 1989. Vol. 315. P. 69.
32. Paul T. Function analitic on half-plane as quantum mechanical states // J. Math. Phys. 1984. Vol. 24. P. 136.
33. Grossman A., Morlet J. Decomposition of hardly functions into square integrable wavelets of constant shape // SIAM J. Math. Anal. 1984. Vol. 15, № 4. P. 273.
34. Sitnikova E.Yu., Hramov A.E., Koronovskii A.A., van Luijtelaar G. Sleep spindles and spike-wave discharges in EEG: Their generic features, similarities and distinctions disclosed with Fourier transform and continuous wavelet analysis // Journal of Neuroscience Methods. 2009. Vol. 180. P. 304.
35. Sitnikova E., van Luijtelaar G. Cortical and thalamic coherence during spike-wave seizures in WAG/Rij rats // Epilepsy Res. 2006. Vol. 71. P. 159.
36. Абдуллаев Н.Т., Дышин О.А., Самедова Х.З. Автоматическая классификация электроэнцефалограмм на основе их вейвлет-пакетной обработки // Биомедицинская радиоэлектроника. 2009. № 6. C. 45.
37. Абдуллаев Н.Т., Дышин О.А., СамедоваХ.З. Вейвлетная очистка электроэнцефалограмм от артефактов с адаптацией к их виду и динамике // Биомедицинская радиоэлектроника. 2009. № 12. C. 34.
38. Божокин С.В., Суворов Н.Б. Вейвлет-анализ переходных процессов электроэнцефалограммы при фотостимуляции // Биомедицинская радиоэлектроника. 2008. № 3. C. 85.
39. Pearson E.S., Neyman J. On the problem of two samples. Joint Statistical Papers, Cambridge University Press, Cambridge, 1967
40. Raiffa H. Decision analysis: Introductory lectures on choices under uncertainty. Addison-Wesley, Reading, 1968.
41. van Luijtelaar G., Hramov A.E., Sitnikova E.Yu., Koronovskii A.A. Spike-wave discharges in WAG/Rij rats are preceeded by delta and theta precursor activity in cortex and thalamus // Clinical Neurophysiology. 2011. Vol. 122. P. 687.
Саратовский государственный Поступила в редакцию 14.03.2011
университет им. Н.Г. Чернышевского После доработки 23.05.2011
Институт высшей нервной деятельности
и нейрофизиологии РАН
WAVELET ANALYSIS OF SLEEP SPINDLES ON EEG AND DEVELOPMENT OF METHOD FOR THEIR AUTOMATIC DIAGNOSTIC
V.V. Grubov, A.A. Ovchinnikov, E.Yu. Sitnikova, A.A. Koronovskii, A.E. Hramov
The detailed wavelet analysis of sleep electric brain activity, obtained from rats with genetic predisposition to absence-epilepsy, has been performed. Characteristic features of time-and-frequency structure of sleep spindles (oscillatory pattern, that serve as electro-encephalographic correlate for slow-wave sleep) have been discovered in long-term
107
electroencephalographic data. Operation has been performed using continuous wavelet transform. Few common wavelet bases have been tested and complex Morlet-wavelet turned out to be the most effective for detection of time-and-frequency features of sleep spindles on EEG. Morlet-wavelet has been used for development of system for automatic diagnostic of sleep spindles on EEG. As a result of analysis two types of sleep spindles, that have the same time dynamics, but different frequency structure, have been discovered. Complex dynamics of main frequency during the sleep spindle has been revealed. The method for automatic diagnostic of sleep spindles, based on computation of wavelet transform energy in two frequency ranges for two types of sleep spindles, has been proposed according to obtained data. The testing of method revealed high accuracy of automatic diagnostic for investigating events on EEG. The method can be used in routine EEG researches, related to detection and classification of different oscillatory patterns.
Keywords: Continuous wavelet transform, EEG, epilepsy, oscillatory pattern, sleep spindles, diagnostic.
Грубое Вадим Валерьевич - родился в пос. Емильчино Житомирской области (Украина). Студент факультета нелинейных процессов Саратовского государственного университета им. Н.Г. Чернышевского. Работает лаборантом в лаборатории № 5 ОФНС НИИ ЕН СГУ
410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н. Г. Чернышевского E-mail: vvgrubov@gmail.com
Овчинников Алексей Александрович - родился в Саратове (1985), окончил Саратовский государственный университет (2008), защитил диссертацию на соискание ученой степени кандидата физико-математических наук (2011, СГУ). Работает в СГУ на кафедре электроники, колебаний и волн в должности ассистента. Область научных интересов: влияние шумов на явление синхронизации в живых системах; диагностика осцилляторных паттернов.
410012 Саратов, ул. Астраханская, 83
Саратовский государственный университет им. Н.Г. Чернышевского E-mail: a.o.ovtchinnikov@gmail.com
Ситникова Евгения Юрьевна - родилась в Волгограде (1974). Окончила ветеринарно-биологический факультет Московской академии ветеринарной медицины и биотехнологии им. К.И.Скрябина (1996), кандидат биологических наук (специальность 03.00.13 - физиология), phD (Social science, Radboud University Nijmegen, Нидерланды, 2008). Старший научный сотрудник Учреждения Российской Академии Наук Института высшей нервной деятельности и нейрофизиологии РАН. Научные интересы - междисциплинарные исследования функций мозга. Автор около 40 научных работ, опубликованных в отечественных и международных журналах.
117485 Москва, ул. Бутлерова, 5А
Институт высшей нервной деятельности и нейрофизиологии РАН E-mail: eu.sitnicova@gmail.com
108
