CC BY
41
Радиофизика
Вестник Нижегородского университета им. Н.И. Лобачевского, 2013, № 1 (1), с. 73-78
УДК 530.182, 57.087
АНАЛИЗ СВЯЗЕЙ МЕЖДУ ОТВЕДЕНИЯМИ ЭЛЕКТРОЭНЦЕФАЛОГРАММ КРЫС ДО И ВО ВРЕМЯ ЭПИЛЕПТИЧЕСКОГО ПРИПАДКА С ПОМОЩЬЮ ПРЕДСКАЗАТЕЛЬНЫХ МОДЕЛЕЙ
© 2013 г. М.В. Сысоева 1, Т.В. Диканев 1, И.В. Сысоев 1, Б.П. Безручко 1,2
'Саратовский госуниверситет им. Н.Г. Чернышевского 2Саратовский филиал Института радиотехники и электроники им. В.А. Котельникова РАН
ivssci@gmail.com
Поступила в редакцию 18.07.2012
Реконструкция работоспособных моделей по экспериментальным временным рядам зачастую требует разработки специфических подходов. В работе с помощью уточненных предсказательных моделей методом нелинейной причинности по Грейнджеру анализируется связанность отведений внутричерепных энцефалограмм крыс - моделей абсанс-эпилепсии. Проведённые расчёты указывают на увеличение связанности за 1.5-2 с до эпилептического разряда. Сценарий изменения уровня связанности в течение припадка различен для разных каналов электроэнцефалограммы.
Ключевые слова: связанность, причинность по Г рейнджеру, эпилепсия, электроэнцефалограмма.
Введение
Определение причинности по Грейнджеру [1] является популярным методом поиска наличия и направления связи между двумя источниками сигналов только по дискретным записям их колебаний (временным рядам). Для анализа причинности по Грейнджеру вначале строится индивидуальная модель, например вида:
Хп+Р=/ ( Хп,Хп-1,...,Хп< Вх _1); ) + 5 п , (1)
учитывающая точки только из одного ряда {хп }Г=1 , влияние на который оценивается. Затем — совместная модель (2), которая учитывает точки из обоих рядов {хп }П=1 и {уп }= :
Хп+Р=Гху (Хп,..,ХпЧ Вх-1)/,Уп,...,Уп-(Ву-1)/ ) + 5п , (2)
где /ху - аппроксимирующая функция, например полином общего вида от Бху = Бх + Бу переменных, Р - дальность прогноза, 5 - остатки модели. Для обеих моделей рассчитываются среднеквадратичные ошибки прогноза: для индивидуальной - 82, для совместной - г2ху . При этом в качестве признака существования связи используется уменьшение ошибки прогноза 82
будущего, удаленного на Р интервалов измерения, значения из первого ряда с помощью модели, построенной по данным только первого сигнала, при введении в модель данных из ряда второго. Количественная мера относительного
улучшения прогноза («причинность по Грейнджеру») имеет вид:
є2 - є2
G = х ху. (3)
єх
В настоящее время эта мера, пришедшая из экономики, применяется для диагностики связанности в различных областях науки, в том числе в нейрофизиологии [2-6] и климатологии [7, 8].
Естественно, что описанный подход не мог остаться незамеченным при решении задач исследования эпилепсии, остающейся одной из актуальнейших проблем физиологии. В работах [9, 10] оценка «причинности по Грейнджеру» применялась для анализа данных записей электроэнцефалограмм (ЭЭГ), содержащих участки (разряды), характерные для больных аб-санс-эпилепсией (рис. 1а). Эта форма заболевания головного мозга проявляется во внезапном очень кратковременном отключении сознания; ею чаще всего страдают дети.
Экспериментальные данные предоставлены нам коллегами из Нидерландов, из университета Наймегена. Они представляют собой набор временных реализаций четырёхканальных внутричерепных ЭЭГ крыс генетической линии WAG/Rij [11], специально выведенных в качестве модели абсанс-эпилепсии. Ряды записаны с частотой дискретизации 1024 Г ц; каждая запись содержит несколько десятков эпилептических припадков различной длины. Каналы соответствуют следующим областям мозга: 1 - затылочная доля коры головного мозга (ОС); 2 -передняя доля коры головного мозга ^С); 3 -ретикулярное ядро таламуса (ЯТЫ); 4 -
вентропостеромедиальное ядро таламуса (УРМ).
В [9] оценка причинности по Грейнджеру показала увеличение связанности во время припадка, а в [10] было обнаружено упрощение динамики, выражавшееся в снижении минимальной достаточной для описания процесса размерности модели непосредственно перед припадком. В [9, 10] использовались линейные модели с единичной дальностью прогноза, а из всех возможных параметров метода подбиралась лишь размерность вектора состояния. Однако опыт применения причинности по Грейнджеру для оценки связанности говорит, что специальный подбор параметров модели может обеспечить большую точность анализа. Поэтому целью нашей работы стало получение более детальной информации об изменениях колебательной активности в фоне до припадка, непосредственно перед его началом и в первое время после его начала с помощью специализированной прогностической модели, использующей неравномерное вложение при конструировании вектора состояния по временному ряду.
Структура предсказательной модели
Реконструкция прогностической модели по временному ряду предполагает выбор способа задания динамической переменой, задания вида оператора эволюции и подбора коэффициентов. Практика показывает, что для успешности метода причинности по Грейнджеру ключевым этапом является выбор формы эмпирической модели. Мы реконструировали модель в виде отображения последования с полиномом общего вида / в качестве аппроксимирующей функции [12], поскольку колебания сложные, нерегулярные, временной ряд короткий (как в смысле числа точек, так и в смысле числа характерных периодов):
Хп+Р=/(хп,хп-1’...,хп-(В-1)1 )+ 5п , (4)
где Р - дальность прогноза, / - лаг, 5п - остатки модели, по смыслу являющиеся шумом. Коэффициенты модели подбирались методом наименьших квадратов.
Величина дальности прогноза выбиралась равной Р = Т / 4 , где Т - характерный период колебаний, поскольку согласно [13-15] этим обеспечивается наибольшая чувствительность метода причинности по Грейнджеру при сохранении малой вероятности ложных выводов. Период определялся по первому наибольшему максимуму на графике автокорреляционной функции (см. рис. 1Ь) либо по спектру мощности (рис. 1с).
Чтобы подобрать порядок полинома V и размерность модели Б, мы руководствовались необходимостью соблюсти баланс между возможностью недообучения модели, когда недостаточная её сложность не позволит в должной мере воспроизвести наблюдаемую динамику, и возможностью переобучения, когда излишне большое число коэффициентов не может быть достаточно достоверно оценено по экспериментальным данным и, как следствие, отражает шумы. Качество модели характеризуется в первую очередь среднеквадратичной ошибкой аппроксимации 8 , равной по построению дисперсии шума 5, однако при условии достаточ-
^ ~ 2 ной точности вычислений 8 всегда является
невозрастающей функцией числа коэффициентов, поэтому с точки зрения минимизации ошибки любое наращивание размера модели оправдано. Баланс между числом коэффициентов и величиной ошибки можно соблюсти воспользовавшись одним из статистических информационных критериев [14, 16], в данной ситуации было решено воспользоваться критерием Шварца [17] - минимумом функции (5):
5 = # • 1п(82 )+С • 1п(#). (5)
Как было показано в работе [18], при использовании стандартной модели вида (4) либо функция 5 вообще не имеет чёткого минимума при моделировании припадка, либо для некоторых записей он достигается при очень большом числе коэффициентов. Там же был предложен способ решения этой проблемы - идея неравномерного вложения (см. также [19]): необходим учёт точки, лежащей на расстоянии порядка характерного периода Т, т. е. в дополнение к точкам, выбранным с лагом /1 =/, следует ввести в вектор состояния одну точку на расстоянии Т, таком, что приближённо выполняется соотношение
Р + /Т=Т. (6)
Для простоты мы ограничились введением только линейного слагаемого в последнем случае, таким образом, в модель (1) добавится всего один дополнительный коэффициент ас+х :
х,„=п х ,х , ,...,х , ч I + аг,,х , + 5,(7)
п+р ./ I п> п-/1> ' п-(вх-1)/1 I С+1 п-/Т Ъп,'~ '
а в совместную модель (2) - два дополнительных слагаемых:
Х„+Р = !ху ^хп,...,хп-(Ох-1)/1 ’Уп’...’Уп-(Ву-1)/1 ^ + (8)
+ аС+1 Хп-/Т + ас+2 уп-/Т + 5п .
В результате для дальности прогноза, равной четверти характерного периода, для всех крыс
0.6
0.4
0.2
я 0
аі -0.2
'ч -0.4
-0.6
-0.8
-1
10
15
20
г, вес (а)
25
30
35
О
Со
0.2 0.4 0.6 0.8
Ґ, 8ЄС
(Ъ)
/,Нг
(с)
Рис. 1. Эпилептический разряд для одного из животных: (а) - временной ряд до, во время и после припадка; (Ь) -автокорреляционная функция, рассчитанная во время припадка (сплошная линия) и до него (пунктир); (с) -спектр мощности, рассчитанный во время припадка (сплошная линия) и до него (пунктир)
были выбраны параметры нелинейной компоненты: порядок полинома у = 3, размерность Бх = 6, лаг /1 = 10 и дополнительный лаг 1Т=Т - Р .
Временные ряды до и во время припадка существенно различаются, характеризуясь высокой степенью нестационарности [20, 21]. Тем не менее детальный анализ автокорреляционной функции (см. рис. 1Ь) показывает, что и до припадка в сигнале имеется регулярная компонента с тем же периодом, что во время, хотя она выражена существенно меньше. Это позволяет до некоторой степени рассчитывать на эффективность построенной в [18] модели для всего временного ряда. Даже если дополнительное слагаемое с коэффициентом ас+1 окажется лишним, полученные в [10] результаты свидетельствуют, что модель всё равно будет приемлемо описывать доразрядную энцефалограмму, поскольку используется достаточно большая размерность Бх = 6.
Поскольку мы рассматриваем вопрос о наличии или отсутствии связи, было решено, что достаточно ограничиться величиной Бу= 1, тем более что некоторые расчёты, проведённые при Dy = 3, дают качественно те же результаты.
Экспериментальные результаты и выводы
Чтобы оценить, как связанность между различными областями мозга меняется во времени, в том числе с началом припадка, величина грейнджеровской причинности вычислялась в скользящем окне длиной 0.5 с со сдвигом между последовательными окнами 0.1 с. Столь малая длина окна была выбрана нами для получения высокого временного разрешения. Следует заметить, что расчёты, проведённые с использованием окон длиной 1 с, показали качественно те же результаты. Для каждого припадка в каждой записи исследовались 10 с непосредственно до и первые 5 с во время него; припадки длиной менее 5 с не рассматривались. Результирующие зависимости G {¡:) усреднялись по всем припадкам для каждого животного отдельно, и среднее выводилось на график. На рис. 2 приведены графики, построенные для всех пар отведений в обе стороны, т. е. для каждой крысы получается 12 графиков, всего - 60. Различными кривыми показаны зависимости G{¡:) для различных животных.
2 —> 3
3 —> 2
и я 2 —>4
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
1 1 1 1 4 і і
_ : : : : Г;'А - : ; : " і \/ і: і ‘ і і ! і і
0.7 1 1
0.65 - і І-
0.6 - : і-
0.55
0.5
0.45
0.4 0 35 .$¿¡¡.>1 .»і!..?.« І: і і
...
ї> \.'\;А і Ш ■
Ху ■'■ :<т
я; :
и $
2 -> 1
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
1 1 1 1 і : і : : А, і і
і і і л.А'- Л» .А »і НІ У$і£Ві г лі * її щ і1 :: і і і *' і ! і і
и я 1 ~>3
¡Л,
Ліі * 1І
. : „і. і 1 л/\ЩГ •
VWУ^У^?W\iД/» г . %
;л.ч,,. ,.^¿4*
\ л'ЬлЛ'\..Л*Г.
..^1%;
уі> лі: л А;
I, я 3-М
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
1 1 1 1 1 : ; ; : 1 1 1
: ; ; : І : ; : : і : ; ; : 1,
; і | д*' ¿4/лй
їШШ [Ч*-г-1 1
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.7 1 1 1 1 і 1 1 0.7
0.65 ■ ' ; ! : ; л|/\ : ; ; . • ■ і} \! Іі ’ - ■ - -А. ■ /\./Цц\ -/лЕр, 2^4^ '/Л К Щ 0.65
0.6 0.55 0.5 0.6 0.55 Й 0.5
0.45 Н 0.45
0.4 і : і- 0.4
0 35 і •' ’ і Iі ' і ! 1 1 0.35
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
І І І І , 4. . : : І І
' : і і і ' і ! і.л ..¥...VI. .....¿ Шг* І і
0 2 4 6 8 10 12 14 і, я 3 -> 1
! ! ! ! К і і і і ! І і
0 2 4 6 8 10 12 14 і, я 4 -> 1
і і і і і - ч і І / "¿5 ї4і''.л-..Xі У ' Пї\/*И У\іт іуя 'у ІЛ’ І 1 •< : /1 .Лі /VI: и і у » у ’у^у І І/и / Ы^/'г: 1 : : V : і і і і і ! І ! ! _
0 2 4 6 8 10 12 14 і, я 4~>3
і і і і і : : : ; 1 : ; : ; і : : Г ; І" ■ і І І і і
Рис. 2. Зависимость улучшения прогноза от времени, построенная в окнах длиной 0.5 с для всех отведений. Зависимости, полученные для различных животных, показаны на одном графике разными кривыми, каждый график соответствует какой-либо паре отведений. Вертикальная сплошная линия обозначает начало припадка, отмеченное физиологами. Расстояние от штриховой линии до сплошной показывает длину скользящего о кна. Его смысл в том, что если изменения G начались в пределах окна (от пунктирной до сплошной линии), они могут быть обусловлены тем, что окно захватывает припадок
Численный анализ показал, что во всех случаях (см. рис. 2) показатель улучшения прогноза увеличивается за 1-2 с перед эпилептическим припадком, что может служить косвенным свидетельством активизации взаимодействия различных областей мозга, причём глобально -каждая область усиливает своё воздействие на все прочие. Этот вывод можно считать достаточно надёжным, поскольку результат повторяется для всех пяти животных без исключения.
Далее существуют два сценария поведения улучшения прогноза во время припадка, общие для всех животных:
1. Показатель G падает с началом припадка до значений того же порядка или даже ниже, что и в фоновом режиме. Это характерно при воздействии на отведение 2 (переднюю долю коры головного мозга). Воздействие на него незадолго до припадка увеличивается, а с началом припадка уменьшается.
2. В течение всего припадка G остаётся на постоянном уровне, до которого улучшение прогноза поднялось перед ним. Это характерно при воздействии на отведение 1 (затылочную долю коры головного мозга) и отведение 3 (ретикулярное ядро таламуса), т. е. воздействие на них во время и незадолго до припадка усиливается.
При воздействии на отведение 4 (вентропо-стеромедиальное ядро таламуса) причинность по Грейнджеру у разных животных ведёт себя по-разному: перед припадком всегда резко возрастает, а затем либо резко падает в самом начале, либо в течение припадка остаётся на постоянном уровне.
В итоге можно сформулировать следующие выводы о том, как полученные в данной работе с помощью модифицированной модели, предложенной в [18], результаты позволяют дополнить известные ранее [9, 10] исследования. Модифицированная нелинейная модель указывает на увеличение взаимодействия между всеми без исключения каналами за 1-1.5 с перед припадком. В работе [10] на тех же промежутках времени было выявлено начальное упрощение динамики, продолжающееся во время припадка. Можно предположить, что увеличение взаимодействия приводит к синхронизации отведений и, как следствие, упрощению динамики, что выражается в снижении размерности модели, необходимой для адекватного описания, и было отмечено в [9].
В работе [9] наблюдается увеличение взаимодействия между всеми отведениями во время припадка, в то время как наши расчёты показывают увеличение воздействия только на отведения 1 и 3, а также на отведение 4 в ряде случаев, т. е. результаты согласуются частично. В то же
время само значение улучшения прогноза G при использовании нелинейных моделей в 2-3 раза больше, т. е. можно предположить, что использование нелинейных полиномов общего вида в качестве аппроксимирующих функций даёт больше возможностей для выявления связанности, что согласуется с результатами [22, 23]. Таким образом, если полученные в [9] с помощью линейных моделей результаты отражают общее увеличение синхронности — линейные авто- и взаимные регрессионные модели связаны с автокорреляционной и взаимной корреляционной функциями, полученные нами с помощью модифицированной нелинейной модели результаты можно интерпретировать как уточняющие более детальную картину взаимодействия между отведениями во время припадка. Надеяться на это позволяет существенная оптимизация структуры модели именно для припадка, проведённая нами в [18], а также хорошая повторяемость результатов на всех пяти животных, хотя этот вопрос ещё нуждается в дополнительном исследовании, в том числе на моделях.
Работа выполнена при поддержке РФФИ, грант 1202-00377, программы «Фундаментальные науки - медицине» и ФЦП «Научные и научно-педагогические кадры инновационной России».
Список литературы
1. Granger C.W.J. Investigating causal relations by econometric models and cross-spectral methods // Eco-nometrica. 1969. V. 37. № 3. P. 424-438.
2. Winterhaldera M., Schelter B., Hesse W. et al. Comparison of linear signal processing techniques to infer directed interactions in multivariate neural systems // Signal Processing. 2005. V. 85. Р. 2137-2160.
3. Roebroeck A., Formisano E., Goebel R. Mapping directed influence over the brain using Granger causality and fMRI // Neuroimage. 2005. V. 25. P. 230-242.
4. Gourevitch B., Le Bouquin-Jeannes R., Faucon G. Linear and nonlinear causality between signals: methods, examples and neurophysiological applications // Biological Cybernetics. 2006. V. 95. P. 349-369.
5. Bernasconi C., Konig P. On the directionality of cortical interactions studied by structural analysis of electrophysiological recordings // Biological Cybernetics. 1998. V. 81. P. 199-210.
6. Hesse W., Moller E., Arnold M., Schack B. The use of time-variant EEG Granger causality for inspecting directed interdependencies of neural assemblies // J. Neurosci. Methods. 2003. V. 124. P. 27-44.
7. Мохов И.И., Смирнов Д.А. Эмпирические оценки воздействия различных факторов на глобальную приповерхностную температуру // Доклады Академии наук. 2009. Т. 426. С. 679-684.
8. Mokhov I. I., Smirnov D. A., Nakonechny P. I. et al. Alternating mutual influence of El. Nino /Southern oscilla-
tion and Indian monsoon // Geophys. Res. Lett. 2011. V. 38. Issue 8. GL 045932.
9. Sitnikova E., Dikanev T., Smirnov D. et al. Granger causality: Cortico-thalamic interdependencies during absence seizures in WAG/Rij rats // Journal of Neuroscience Methods. 2008. V. 170. P. 245-254.
10. Филина Е.В. Динамика локальных потенциалов мозга при абсанс-эпилепсии: эмпирическое моделирование // Изв. вузов. ПНД. 2011. Т. 19. № 4. C. 109-123.
11. Coenen A.M.L. and van Luijtelaar E.L.J.M. Genetic animal models for absence epilepsy: A review of the WAG/Rij strain of rats // Behavior Genetics. 2003. V. 33(6). P. 635-655.
12. Chaos and its Reconstruction / Edited by G. Gouesbet and S. Meunier-Guttin-Cluzel. Nova Science Publishers, 2003. P. 215-243.
13. Сысоева М.В., Диканев Т.В. Подбор дальности прогноза, оптимальной для оценки связи по Г рейнджеру // Материалы VI конференции молодых учёных «Наноэлектроника, нанофотоника и нелинейная физика». Саратов, 2011. С. 162-163.
14. Безручко Б.П., Смирнов Д.А. Математическое моделирование и хаотические временные ряды. Саратов: ГосУНЦ «Колледж», 2005. 320 с.
15. Сысоева М.В., Диканев Т.В., Сысоев И.В. Выбор временных масштабов при построении эмпирической модели // Изв. вузов. ПНД. 2012. Т. 20. № 2.
16. Nonlinear phenomena research perspectives / Ed. C.W. Wang. New York: Nova Science Publishers, 2007. P. 7-53.
17. Schwarz G. Estimating the dimension of a model // The Annals of Statistics. 1978. V. 6. № 2. Р. 461-464.
18. Сысоева М.В., Сысоев И.В. Математическое моделирование динамики энцефалограммы во время эпилептического припадка // Письма в ЖТФ. 2012. Т. 38. Вып. 3. С. 103-110.
19. Judd K., Mees A. Embedding as a modeling problem // Physica D. 1998. 120. Р. 273-286.
20. Dikanev Т., D. Smirnov, Wennberg R., et al. EEG nonstationarity during intracranially recorded seizures: statistical and dynamical analysis // Clinical Neurophysiology. 2005. V. 116. P. 1796-1807.
21. Глушкова М.В., Диканев Т.В. Отличия временного ряда внутричерепной ЭЭГ крыс в начале и конце эпилептического разряда // Материалы XII Всероссийской школы-семинара «Волновые явления в неоднородных средах». Звенигород, 2009. С. 19-21.
22. Yonghong Chen, Govindan Rangarajan, Jianfeng Feng, Mingzhou Ding. Analyzing multiple nonlinear time series with extended granger causality // Physics Letters A. 2004. V. 324. Issue 1. P. 26-35.
23. Сысоев И.В., Караваев А.С., Наконечный П.И. Роль нелинейности модели в диагностике связей при патологическом треморе методом грейнджеровской причинности // Изв. вузов. ПНД. 2010. Т. 18. № 4. С. 81-88.
ANALYSIS OF COUPLING BETWEEN RAT ELECTROENCEPHALOGRAM CHANNELS BEFORE AND DURING EPILEPTIC SEIZURE USING PREDICTIVE MODELS
M. V. Sysoeva, T. V. Dikanev, I. V. Sysoev, B.P. Bezruchko
The reconstruction of efficient models from experimental time series often requires special approaches. The coupling between rat intracranial electroencephalogram channels (as models of absence epilepsy) is analyzed using improved predictive models by means of nonlinear Granger causality. The calculations have shown the coupling growth 1.5-2 seconds before an epileptic seizure. The scenario of coupling level changes during the seizure is not the same for different electroencephalogram channels.
Keywords: coupling, Granger causality, epilepsy, electroencephalogram.
