CC BY
13
Прикладные задачи
^^^^^^^^^^»нелинейной теории колебаний и вслн
УДК 530.182, 51-73
Оптимизация параметров метода причинности по Грейнджеру для исследования лимбической эпилепсии
М. В. Сысоева1, Т. М. Медведева2
Саратовский государственный технический университет имени Гагарина Ю.А.
Россия, 410054 Саратов, Политехническая, 77 2 Институт высшей нервной деятельности и нейрофизиологии РАН Россия, 117485 Москва, Бутлерова, 5А E-mail: bobrichek@mail.ru, Golova93tanya@gmail.com Автор для переписки Сысоева Марина Вячеславовна, bobrichek@mail.ru Поступила в редакцию 29.03.2018; после доработки 27.05.2018
Цель. Выявить зависимость результатов анализа связанности между отделами лимбической системы мозга, полученных с применением причинности по Грейнджеру, от выбранных временных масштабов эмпирических математических моделей, построенных по временным рядам внутричерепных электроэнцефалограмм, записанных во время лим-бических эпилептических разрядов.
Методы. Используется сочетание методов анализа связанности по экспериментальным рядам и подходов моделирования из первых принципов, основанных на воспроизведении основных временных и частотных свойств экспериментальных сигналов. Такая комбинация для исследования связанности между отделами мозга по внутричерепным электроэнцефалограммам применяется впервые, а именно, в данной работе - для анализа связанности при лимбической эпилепсии, вызванной у крыс линии WAG/Rij введением агониста эндоканнабиноидных рецепторов.
Результаты. В ансамблях из четырёх связанных осцилляторов ван дер Поля с потенциалом Тоды и жёстким возбуждением, Хиндмарш-Розе и ФитцХью-Нагумо получены режимы, воспроизводящие ряд спектральных и амплитудных характеристик сигналов локальных потенциалов при лимбических разрядах. Выбраны оптимальные с точки зрения сочетания чувствительности и специфичности параметры метода причинности по Грейн-джеру. Используя эти параметры, по экспериментальным данным крыс линии WAG/Rij выявлено значимое усиление воздействия со стороны затылочной коры на гиппокамп во время лимбических разрядов примерно за 2 с до начала приступа, ослабевающее до фонового уровня ровно в момент конца разряда.
Обсуждение. Надёжность выводов о связанности является ключевой проблемой применения метода причинности по Грейнджеру к анализу экспериментальных данных. Повышение чувствительности и специфичности метода возможно различными способами, включая увеличение объёма экспериментальных данных и адаптацию параметров метода к спектральным свойствам сигнала, но ни один из них не решает проблему окончательно. Существенно продвинуться в этом направлении, на наш взгляд, позволяет предлагаемый подход, основанный на построении ансамблей осцилляторов, генерирующих сигналы, качественно схожие с экспериментальными.
Ключевые слова: временные ряды, реконструкция прогностических моделей, нестационарные сигналы, анализ связанности, причинность по Грейнджеру, электроэнцефалограмма, лимбическая эпилепсия.
https://doi.org/ 10.18500/0869-6632-2018-26-5-39-62
Образец цитирования: Сысоева М.В., Медведева Т.М. Оптимизация параметров метода причинности по Грейнджеру для исследования лимбической эпилепсии // Известия вузов. Прикладная нелинейная динамика. 2018. Т. 26, № 5. С. 39-62. https://doi.org/ 10.18500/0869-6632-2018-26-5-39-62
Optimization of Granger causation method parameters for the study of limbic epilepsy
M. V Sysoeva1, T.M. Medvedeva2
1Yuri Gagarin State Technical University of Saratov 77, Politechnicheskaya str., 410054 Saratov, Russia 2 Institute of Higher Nervous Activity and Neurophysiology RAS 5A, Butlerova str., 117485 Moscow, Russia E-mail: bobrichek@mail.ru, Golova93tanya@gmail.com Correspondence should be addressed to Sysoeva Marina V., bobrichek@mail.ru Received 29.03.2018; revised 27.05.2018
Purpose. The aim is to reveal the dependence of Granger causality results on chosen time scales of constructed empirical models in application to the task of investigation of evolution of coupling between brain areas during limbic seizures.
Methods. We use combination of methods for coupling analysis of the experimental time series and approaches to modeling from the first principles, which reproduce the main time and frequency properties of the experimental signals. Such a combination use is novel for investigation of the connectivity between the brain areas from intracranial electroencephalogram In this paper, it is used for connectivity analysis in limbic epilepsy provoked in WAG/Rij rats by the introduction of endocannabinoid receptor agonist.
Results. In ensembles of four coupled van der Pol oscillators with the Toda potential and hard excitation, Hindmarsh-Rose systems and FitzHugh-Nagumo systems we found regimes reproducing spectral and amplitude characteristics of the series of local potentials at the limbic seizures. Optimal method parameters were selected to target both sensitivity and specificity of Granger causality. Using these parameters, a significant increase in coupling was detected in the experimental data of WAG/Rij rats from the occipital cortex to the hippocampus during limbic seizures approximately 2 s before the seizure onset. The coupling return to background level immediately after the seizure termination.
Discussion. Reliability of coupling detection procedure outcomes is a key problem in applying the Granger causality method to experimental data. Increasing the method sensitivity and method specificity is possible in various ways, including increasing experimental data amount and adapting method parameters to the signal spectral properties, but none of these approaches solves the problem completely. In our opinion, the proposed approach, based on the construction of oscillator ensembles generating signals qualitatively similar to experimental ones, allows us to make significant progress in this direction.
Key words: time series, predictive models reconstruction, nonstationary signals, connectivity analysis, Granger causality, electroencephalogram, limbic epilepsy.
https://doi.org/ 10.18500/0869-6632-2018-26-5-39-62
References: Sysoeva M.V., Medvedeva T.M. Optimization of Granger causation method parameters for the study of limbic epilepsy. Izvestiya VUZ, Applied Nonlinear Dynamics, 2018, vol. 26, no. 5, pp. 39-62. https://doi.org/ 10.18500/0869-6632-2018-26-5-39-62
Введение
Данная работа является частью большого исследования, посвященного комбинированию методов математического моделирования из первых принципов (прямое моделирование) и метода построения модели путём решения обратной задачи (обратное моделирование) [1] для лучшего понимания того, как работает мозг в норме и при патологиях. Ставится цель выявления зависимости результатов анализа связанности отдельных областей мозга методом причинности по Грейнджеру [2] от выбранных временных параметров предсказательных математических моделей.
В качестве экспериментальных данных используются записи внутричерепных электроэнцефалограмм (ЭЭГ) крыс, содержащие лимбические разряды, вызванные введением агониста эндоканнабиноидных рецепторов [3]. Так как реальный эксперимент является существенно более сложным, дорогостоящим и непредсказуемым этапом работы, чем создание математической модели, то вся работа и отталкивается, в первую очередь, от реальных данных, полученных на экспериментальных животных.
Поскольку взаимодействия между отделами мозга играют главную роль в механизмах протекания эпилептических разрядов [4-9], возникает потребность в использовании современных методов детектирования наличия и направления связей. Метод причинности по Грейнджеру как раз позволяет установить наличие и направленность влияний систем друг на друга [10]. Для этого необходимо, чтобы системы порождали некий процесс, который можно зарегистрировать. Для наших реальных данных этим процессом является электрическая активность мозга, которая регистрируется в виде электроэнцефалограммы. С точки зрения нелинейной динамики электроэнцефалограмма есть не что иное, как временной ряд. Основная же идея метода причинности по Грейнджеру заключается в следующем: если прошлые значения одного временного ряда, измеренного от первой системы, помогают точнее предсказывать будущие значения другого ряда, полученного от второй системы, то считается, что первая система влияет на вторую. Чтобы провести такой анализ, для имеющихся временных рядов нужно построить предсказательные математические модели. В своей оригинальной работе Грейнджер использовал только линейные авторегрессионные модели. В настоящее время успешно применяются более сложные нелинейные модели [11,12], в том числе и в задачах нейрофизиологии [13-20].
Выбор подходящих параметров модели (параметризация) очень важен для успеха метода: даже в линейной грейнджеровской причинности выбор размерности модели (число точек в прошлом, которые используются для предсказания будущего состояния) имеет большое влияние на предсказательную способность [21-24]. Ошибочная параметризация может стать причиной ложных результатов: будут детектированы связи, которых в реальности нет (ошибка первого рода, плохая специфичность метода), или не будут обнаружены реально существующие связи (ошибка второго рода, плохая чувствительность метода) [11,25,26].
Далеко не все параметры можно подобрать с помощью объективных критериев. Обратная задача и вовсе в общем случае некорректна (в том числе из-за ограниченности объёма экспериментальных и априорных данных об изучаемой системе) и имеет множество решений. Ко всему прочему большинство известных статистических критериев применимы только для построения прогностических моделей одной системы, а не для анализа связанности нескольких систем, хотя при реконструкции систем с запаздыванием имеется удачный опыт, когда критерий, разработанный для
автономных осцилляторов [27], оказался применим и для их сетей [28]. Как было показано в [24], параметры, оптимальные для индивидуальной модели, не всегда являются оптимальными для метода грейнджеровской причинности. Поэтому метод нуждается в дополнительном тестировании, без чего выводы о наличии, направлении и изменении связанности оказываются ненадёжны.
Ещё одной проблемой для корректной работы разрешённого во времени метода причинности по Грейнджеру являются быстрые переходные процессы. Как было показано в [29], метод имеет артефакт (феномен «уши») при захвате скользящим окном быстрых переходных процессов (в случае экспериментальных данных такой эффект возникает при переходах от преиктальной к иктальной фазе и от иктальной к постиктальной фазе). Данный артефакт выражается в большинстве случаев в резком сильном возрастании (иногда в падении) причинности по Грейнджеру в начале и в конце «разряда» в пределах ширины скользящего окна. Такое резкое изменение есть следствие переходного процесса и объясняется тем, что построить хорошую индивидуальную модель во время переходного процесса очень сложно из-за нестационарности, в данном случае вызванной быстрым изменением силы связи или изменением собственных параметров эталонного осциллятора. Сложность описания неавтономных и переходных процессов автономными моделями уже неоднократно отмечалась в литературе [30-32]. Таким образом, предсказательная сила собственной модели падает, в то время как вклад добавки, учитывающей влияние второго ряда, основывается на измеренных значениях второго ряда и не зависит от смены индивидуальных параметров. Поэтому даже при падении коэффициента связи, относительный вклад добавки в прогноз может существенно вырасти.
Тестировать метод причинности по Грейнджеру непосредственно на экспериментальных данных невозможно, поскольку для этого нужно знать реальную архитектуру связей. Имеется в виду, в первую очередь, связанность с точки зрения физиологии и теории информации, поскольку связанность с точки зрения морфологии уже достаточно хорошо изучена. Для такого тестирования в данной работе предлагается прибегать к эталонным системам, которые должны при достаточной простоте воспроизводить основные свойства экспериментальных сигналов, по крайней мере, те, которые существенны и учитываются при поиске связанности. В данном исследовании для нас важно соблюсти временные масштабы, характерные для экспериментальных данных (частоту дискретизации и основную частоту колебаний), необходимо учесть факт наличия высших гармоник в спектре, а также воспроизвести переход от низкоамплитудной шумоподобной фоновой активности к высокоамплитудной регулярной эпилептиформной активности.
1. Методика
1.1. Экспериментальные данные. Эксперименты были выполнены нидерландскими коллегами под руководством доктора C.M. van Rijn в Donders Centre for Cognition, Radboud University Nijmegen в Нидерландах, были одобрены местным комитетом по этике (RUDEC 2006-064) и соответствовали требованиям европейского соглашения (European Communities Council Directive, 86/609/EEC).
Данные записывались с поверхности коры головного мозга в затылочной доле [AP-7;L 6] и из лимбической системы в гиппокампе [AP-3.5;L+2, глубина —3] интракраниально. Референсный электрод размещался над мозжечком. ЭЭГ была отфильтрована в полосе пропускания 0.1-10 Гц, оцифрована с частотой дискретизации 512 Гц 16-разрядным АЦП со встроенным усилителем и с полудиапазоном 10 B.
10 5 0 -5 -10
i
10 5 0 -5 -10
Hp (Time series)
0 10 20 30 40 50 60 70 80 90 ОС (Time series)
Щ
|Ч|||||||(||||1Р1|1Р*|И : P"
30 £ 25 20
& g 15
Oi £
Hp (Spectrogram)
. 11, hi If' 1 /п 1 ■
0 10 20 30 40 50 60 70 80 90
ОС (Spectrogram)
10 20 30 40 50 60 70 80 90
t, s
10 20 30 40 50 60 70 80 90
t, s
Рис. 1. Пример записи ЭЭГ и спектрограмм лимбического разряда из гиппокампа (Hp) и затылочной доли коры (ОС), исследованные интервалы ЭЭГ до, во время и после разряда (начало и конец разряда отмечены чёрными вертикальными линиями)
Fig. 1. Example of EEG recordings and spectrograms for the limbic seizure from the hippocampus (Hp) and occipital cortex (ОС). Investigated EEG intervals before, during and after the seizure are shown. The seizure onset and the seizure termination are marked by black vertical lines
Животным вводилось 12 мг/кг агониста эндоканнабиноидных рецепторов \УВД55,212-2. Через три часа после введения препарата развивались эпилептические лимбические судороги. Всего было зарегистрировано 13 эпилептических разрядов.
На спектрограммах, построенных по экспериментальным данным, наблюдается следующая типичная динамика (рис. 1). Чёрные вертикальные линии показывают начало и конец разряда (эти маркеры были получены совместно с физиологами). Началом разряда считается момент увеличения амплитуды сигнала на ЭЭГ в два раза по сравнению с фоном, концом разряда - момент уменьшения амплитуды, соответственно. Дополнительная коррекция разметки проводилась с помощью спектрограмм: во время разряда на ней хорошо видна основная гармоника и другие высокие частоты, которых нет в фоне. До и после лимбического разряда в обоих отведениях нет какого-либо основного ритма. Во время разряда в гиппокампе есть слабый низкоамплитудный сигнал с основной частотой около 4 Гц, в затылочной коре - многочастотный высокоамплитудный сигнал с основной частотой также в районе 4 Гц.
1.2. Данные численного эксперимента. Тестовые осцилляторы были сгенерированы примерно с той же основной частотой колебаний, записаны с той же частотой дискретизации, что и экспериментальные сигналы. Чтобы воспроизвести основные возможные ситуации (однонаправленную связь, двунаправленную связь и отсутствие взаимодействия), для каждого типа систем было сделано по 4 осциллятора, связанных следующим образом: первый воздействует на второй, а третий и четвёртым взаимно влияют друг на друга. Тестовые системы демонстрируют поведение, характерное для реальных сигналов: вначале идёт режим с небольшой амплитудой и шумоподобным спектром («фон»), затем идёт высокоамплитудный сигнал с хорошо выраженной основной частотой («разряд»), потом опять низкоамплитудный без основной частоты колебаний режим (снова «фон»). Переход из одного режима
в другой осуществляется следующими способами: в паре первый-второй осциллятор увеличивается один из параметров первой системы и коэффициент связи, в паре третий-четвёртый осциллятор увеличивается только коэффициент связи. Для каждой системы было сгенерировано по тринадцать «разрядов», чтобы при статистическом анализе результатов воспроизвести экспериментальную ситуацию.
В качестве тестовых систем были выбраны следующие.
1.2.1. Система ФитцХью-Нагумо
Хг = Хг(йг - Хг)(Хг - 1) - у + 1а,г + кх3, (1)
Уг = ЬгХг - угуг +
Коэффициент связи к менялся от 0.01 в фоне до 0.4 во время разряда, во время разряда 1а,\ возрастало до значения 1.08, прочие параметры - в таблице 1. Временной ряд и спектрограмма одного из «разрядов» первого осциллятора системы приведены на рис. 2, a.
Таблица 1 (Table 1)
Параметры связанных систем ФитцХью-Нагумо Parameters of coupled FitzHugh-Nagumo systems
i 1 2 3 4
j 1 4 3
ai 0.8 0.8 0.8 0.8
bi 0.15 0.17 0.15 0.17
Y i 0.06 0.068 0.06 0.068
ia,i 0.85 0.86 0.85 0.86
1.2.2. Система Хиндмарш-Розе
Xi = Vi - aix3 + bix2 - zi + Ia,i + kx2 + ,
Vi = ci- dix2 - Vi, (2)
zi = ri (si(xi XR,i) zi).
Коэффициент связи k менялся от 0.01 в фоне до 0.65 во время разряда, во время разряда Ia,i возрастало до значения 4, прочие параметры приведены в таблице 2. Временной ряд и спектрограмма одного из «разрядов» первого осциллятора системы приведены на рис. 2, b.
Таблица 2 (Table 2)
Параметры связанных систем Хиндмарш-Розе Parameters of coupled Hindmarsh-Rose systems
i 1 2 3 4
j 1 4 3
ai 0.55 0.55 0.52 0.5
bi 4.8 4.5 4.5 4.2
Ci 1 1 1 1
di 7 7 7 7
si 6 6 6 6
Ti 10-3 10-3 10-3 10-3
XR,i -1.6 -1.6 -1.6 -1.6
ia,i -2.9 -2.86 -2.8 -2.86
1.2.3. Обобщённый осциллятор ван дер Поля-Тоды
Хг - (г'г + кх] - xf) ±i + w2 (1 - ехр(ж»)) = ^¿(i), (3)
где г - номер текущего осциллятора, j - номер воздействующего осциллятора, коэффициент связи к менялся от 0.01 в фоне до 0.65 во время разряда, во время разряда Г\ возрастало до значения 1.08, '£,;(/) нормальный белый шум. Параметры для всех четырёх систем приведены в таблице 3. Временной ряд и спектрограмма одного из «разрядов» первого осциллятора системы приведены на рис. 2, с.
Таблица 3 (Table 3)
Параметры связанных осцилляторов ван дер Поля-Тоды Parameters of coupled van der Pol-Toda oscillators
г 1 2 3 4
j 1 4 3
Гг -0.08 -0.14 -0.12 -0.11
Wi 0.4 0.4 0.4 0.4
&20 g 15 g. 10 щ 5 § 0
Hindmarsh-Rose (Time series)
x 25
Hindmarsh-Rose (Spectrogram)
Рис. 2. Временные ряды и спектрограммы тестовых систем: а - ФитцХью-Нагумо, Ъ - Хиндмарш-Розе, с - ван дер Поля с жёстким возбуждением и потенциалом Тоды. На рисунке приведены исследованные интервалы ЭЭГ до, во время и после «разряда» (начало и конец «разряда» отмечены чёрными вертикальными линиями)
Fig. 2. Time series and spectrograms of test sistems: a - FitzHugh-Nagumo, b - Hindmarsh-Rose, с - van der Pol with hard excitation and Toda potential. Investigated EEG intervals before, during and after «seizure» are shown. The «seizure» onset and the «seizure» termination are marked by black vertical lines
1.3. Анализ связанности методом причинности по Грейнджеру. При анализе связанности методом причинности по Грейнджеру вначале строится индивидуальная модель, учитывающая данные только из одного ряда {хпвлияние на который оценивается по формуле
п+т
=f (
хп Хп — 1х
1-(Пе-1)1х) + аг., + 1Хп-1тх
(4)
где хП+т - предсказанное значение, соответствующее измеренному значению хп+т; f - полином общего вида порядка Р от Д, переменных [11]; Хп = (хп,хп-1х,..., хп-(рв-1)1х) - вектор состояния, полученный методом задержек [33]. Метод задержек является классическим подходом для получения высокоразмерного вектора состояния {хп}п=1 из скалярного временного ряда {хп}«=1 путем сдвига назад во времени (Дв — 1 )1Х раз на временную задержку 1Х. Здесь: 1тх - добавочный лаг, который учитывает значение экспериментальных данных, задержанное от предсказываемого на величину характерного периода Т; т - дальность прогноза, то есть расстояние между последней точкой в векторе состояния и предсказываемой (на сколько точек вперед предсказываем). Коэффициенты модели подбирались методом наименьших квадратов [34].
Затем строится совместная модель, которая учитывает данные из двух рядов -ряда предположительно ведомой системы {хп}«=1 и ряда предположительно ведущей системы {уп}«=1,
п+т
= 9 (хп,хп-1х >... ,хп-(д3-1)1х ,Уп,Уп-1у,... >Уп-(£а-1)гу) +
+ аХз +1хп-1тх + +2Уп-1ту
(5)
где хп+т - прогнозируемое значение для экспериментального значения хп+т, полученное с использованием значений двух временных рядов {хп}«=1 и {уп}«=1; Д- -размерность совместной модели, равная сумме размерности индивидуальной модели Д3 и добавочной размерности Да. Тогда 9 - полином общего вида порядка Р от Д- = Дв + Да переменных, 1У - лаг для вектора состояния из ряда {уп}«=1.
На основе ошибок прогноза моделей (4) и (5) составляется основная мера связанности Р1 - показатель улучшения прогноза.
1
N £
М'о2 ^ ^хп х'п)
Х n=(Ds-1)lx+ т
1
£; =
Р1
1х
N £
II №
Пр _ Пр
п
N02 ^ У п "Ьп
Х п=(ша.х^^а)-1)1х+т £2
1---
1 £2 '
(6)
(7)
(8)
где хп - измеренное значение; хп - значение, предсказанное индивидуальной моделью; хп - значение, предсказанное совместной моделью; оХ - дисперсия временного ряда (хп}„=1; N' - эффективная длина временного ряда {х п }n=1, вычисляемая как N' = N -((Д, — 1 )1х + т).
2
2
£
в
Если е2 = е2, то есть Р1 = 0, значит сигнал {уп}^=1 никак не помогает предсказать динамику |хга}^=1 и делается вывод, что вторая система не влияет на первую. Если же е2 > 0 и е2 ^ 0, значит Р1 ^ 1, временной ряд {уп}^=1 существенно помогает в прогнозе ряда {хга}^=1 и делается вывод, что вторая система влияет на первую. Отсюда, значение улучшения прогноза Р1 лежит в диапазоне [0; 1]. Теоретически ситуации Р1 = 0 и Р1 = 1 достижимы, но для этого структура моделей (4), (5) должна очень хорошо описывать объект, как это показано в [24], иначе при наличии связи и при её отсутствии будет получено значение 0 < Р1 < 1, что и происходит на практике. Это может быть обусловлено конечностью разрешения по времени [35], неверным выбором параметров метода, конечной длиной ряда, шумами и иными факторами.
Необходимо учитывать, что абсолютное значение величины Р1 несёт мало информации о степени связанности подсистем, как это показано в [36]. Однако увеличение или уменьшение Р1 при использовании одной и той же модели на всём протяжении исследуемого ряда имеет смысл. Таким образом, можно отслеживать изменения силы связей, если анализ связанности проводить в скользящем окне, как это предложено в [37].
1.4. Параметризация метода. Хотя в [23] было показано, что оптимальной дальностью прогноза является четверть характерного периода колебаний, в дальнейших работах [38] высказывалось предположение, что некоторые другие значения могут давать лучшие результаты для конкретных данных. Поэтому в данном исследовании рассматривалось несколько значений дальности прогноза т = 1, Т/16, Т/12, Т/8, Т/4, из которых в дальнейшем, опираясь на результаты, полученные на тестовых системах, выбирали оптимальное значение для имеющихся у нас реальных данных.
Так как для реальных систем очень важно временное разрешение полученных результатов, то влияние реальных и тестовых систем оценивалось переменным во времени адаптированным методом причинности по Грейнджеру [37]. Все расчёты производились в скользящем окне. Ширина скользящего окна и> выбиралась равной 1 с (4 характерных периода колебаний) и 2 с (8 характерных периодов колебаний).
Для каждой тестовой системы и реальных данных по критерию Шварца [39] были подобраны оптимальные параметры прогностических математических моделей (порядок степенного полинома, размерность, основной и добавочный лаг модели) для «фона» и «разряда»
5 = N Мей + г^, (9)
где г - число коэффициентов модели. В классическом случае эта величина вычисляется для индивидуальной модели, как = (Р + В3)\/(РШ81). В нашем случае количество коэффициентов индивидуальной модели увеличивается ещё на один линейный член г = гз + 1.
Так как в данной работе, в первую очередь, ставится задача верного определения направления связи между двумя системами, а не выбора наилучшей индивидуальной прогностической модели, то в критерии Шварца есть смысл закладывать количество коэффициентов совместной, а не индивидуальной модели, а значит, необходимо учесть добавочную размерность ^ и ещё один линейный член. Тогда
количество коэффициентов будет равно г? = (Р + + Ба)!/(Р !(Д, + ^а)1) и ещё два линейных члена г = г? + 2. Берётся только тот набор параметров, который будет удовлетворять условию г ^ д/Ж'.
1.5. Оценка прогностической силы использованных эмпирических моделей. В работе [24] было показано, что хорошая чувствительность и специфичность метода причинности по Грейнджеру могут быть получены при использовании моделей, далёких от идеальных (дающих существенно большую ошибку прогноза). С другой стороны, в [11] показано, что недостаточно хорошие модели (дающие слишком большую ошибку) заметно снижают чувствительность метода. Поэтому важно оценить, насколько ошибка прогноза совместной модели близка к предельно возможной для выбранной дальности прогноза. Минимально возможная ошибка е^ может быть рассчитана по формуле
п2
е^ш = е^1Т (10)
сией шумов оПо1зе измерения, незнания и динамического шума, как это подробно
и определяется старшим ляпуновским показателем системы Xi и суммарной диспер-J2
ПС
изложено в [1].
Поскольку шумы незнания для экспериментальных данных не доступны измерению и даже оценке, полученная величина будет явно занижена. Тем не менее, приближённо можно сказать, что если е^щ того же порядка, что и е2 или на 1-2 порядка меньше, то модель явно хорошая и чувствительность должна быть достаточно высокая. Если е^т на несколько порядков меньше е2, то модель нуждается в улучшении структуры и чувствительность может быть недостаточна. Если же е^щ ~ е2, то модель, скорее всего, переобучена и описывает локальные шумы, а не лежащие в основе изучаемых явлений закономерности.
Чтобы сопоставить полученную ошибку прогноза совместной модели е2 с минимально возможной ошибкой е^щ, рассчитываемой по формуле (10) и таким образом проверить адекватность построенных моделей, необходимо было оценить старший ляпуновский показатель Xi. Для этого был использован метод, предложенный в [40] и зарекомендовавший себя как наиболее устойчивый для малых длин рядов. Расчёт производился по временному ряду, что было необходимо, поскольку во всех модельных системах фигурировал динамический шум, а для экспериментальных данных оператор эволюции просто неизвестен.
Для модельных систем при оценке o2oise использовались априорные знания о включённых в уравнения динамических шумах. Для экспериментальных данных можно приблизительно оценить только измерительный шум, составляющий порядка 4 младших разрядов АЦП. Такая оценка явно существенно занижена и приводит к очень малому отношению етп/е2. Но в действительности, наверняка её чувствительность и специфичность должны быть лучше из-за недооценки шумов при расчёте минимальной ошибки прогноза.
Далее для выбранных значений т = 1, T/16, T/12, T/8, T/4 анализировалось отношение е^щ^. Если это отношение больше единицы, то модель переобучена, при её использовании у метода причинности по Грейнджеру будет плохая
специфичность. Если это отношение очень маленькое (порядка 10-4 и меньше), то модель обладает плохой предсказательной способностью, а метод причинности по Грейнджеру будет иметь плохую чувствительность.
1.6. Статистический анализ полученных результатов. Для тестовых и реальных данных процедура статистического анализа полученных значений улучшения прогноза выглядела одинаково. Наличие 13 реализаций в каждом численном эксперименте (для реальных и тестовых данных) дало возможность произвести статистический анализ значимого отличия полученных значений Р1 (£) от фонового уровня. Фоновый уровень Р^ вычислялся как среднее по всем реализациям по первым 7 секундам. Значения Р1 для каждого момента времени £ рассматривались как выборка. С помощью одновыборочного двустороннего ^теста Стьюдента [41] оценивались значимые отличия среднего этой выборки Р/теап от фонового уровня Р/ь§ на уровне значимости 5%.
На графиках зависимости улучшения прогноза Р1теап от времени £ с оценкой значимости полученных результатов черные точки соответствуют значению Р1теап значимо большему значения фонового уровня Р^ на уровне значимости 0.05. Над каждым графиком написано влияние какой системы на какую проверяется у ^ х.
2. Результаты и обсуждение
2.1. Результаты параметризации. По критерию Шварца, учитывая вышеописанные ограничения, для всех тестовых примеров подобрались одинаковые параметры: Р = 3, Д, = 2, 1 = Т/12, 1т = Т — т. Для реальных данных потребовалась меньшая нелинейность, гораздо важнее оказалась большая размерность: Р = 2, Д, = 4, 1 = Т/12, 1т = Т — т.
Для полученных эмпирических моделей была оценена их прогностическая сила (отношение £¡^/£2), рассчитанная описанным в разделе 1.5 способом. Старший ляпуновский показатель вычислялся для всех «разрядов» по 10-секундному интервалу (около 40 характерных периодов колебаний), вычислялось среднее и среднеквадратичное отклонение полученной выборки. На основе полученных результатов составлена таблица 4.
При единичной дальности прогноза т = 1 для осциллятора ван дер Поля-Тоды 41п/£2 ~ 11.74, для системы ФитцХью-Нагумо £¡¡^/£2 ~ 12.17. Это указывает на то, что модель переобучена, соответственно, можно ожидать плохую специфичность метода причинности по Грейнджеру при таких параметрах модели. Для остальных значений дальности прогноза отношение £¡¡^/£2 порядка « 10-2. При такой модели в методе причинности по Грейнджеру будет наблюдаться оптимальное соотношение чувствительности и специфичности. Для системы Хиндмарш-Розе при всех использованных дальностях прогноза ^ш/^ порядка « 10-4. Это говорит о плохой чувствительности метода грейнджеровской причинности при использовании такой прогностической модели. Скорее всего, в случае системы Хиндмарш-Розе использование короткого скользящего окна не позволяет сделать модель высокоразмерной, поэтому и чувствительность метода недостаточная.
Для экспериментальных данных старший ляпуновский показатель получился равен « 12, а отношение £¡^/£2, как и ожидалось, - заниженным. Можно предпо-
Таблица 4 (Table 4)
Прогностическая сила эмпирических моделей (отношение е^ш/е?)* Predictive power of empirical models (relation е^ш/е?)**
Система ФитцХью-Нагумо Хиндмарш-Розе Ван дер Поля-Тоды Гиппокамп Затылочная кора
Xi, c-1 3.8±0.65 6.2±0.42 8.9±0.59 12.2±0.79 12.5± 1.15
o2 b2 °noise> B 0.02 • At 0.25 • At 4 • At 2.38 • 10-5 2.38 • 10-5
„X, B2 0.13 3.11 1.70 0.28 8.14
T, c e2- /e2
At 0.15 15.9-10-5 25.6-10-4 8.7-10-5 8.740-5
2.9-10-6 =11.74(*) 0.17 =93.4-10-6(o) 2.1-10-4 =12.17(*) 2.540-3 =0.03 2.9-10-3 =0.03
TAt/16 0.16 17.3-10-5 28.9-10-4 1.040-4 1.040-4
2.9-10-3 =0.12 0.43 =40.2-10-6(o) 8.2-10-2 =0.04 0.06 =1.740-3 0.21 =0.540-3
TAt/12 0.17 19.M0-5 33.2-10-4 1.240-4 1.240-4
14.2-10-3 =0.03 0.56 =34.0-10-6(o) 15.340-2 =0.02 0.11 =1.M0-3 0.42 =0.340-3
TAt/8 0.18 20.5-10-5 36.840-4 1.440-4 1.440-4
24.M0-3 =0.02 0.61 =33.6-10-6(o) 20.640-2 =0.02 0.16 =0.9-10-3 0.46 =0.3-10-3
TAt/4 0.19 23.M0-5 43.8-10-4 1.840-4 1.840-4
35.1-10-3 =0.01 0.72 =32.M0-6(o) 25.740-2 =0.02 0.25 =0.740-3 0.47 =0.4-10-3
* Если рядом с полученным числом стоит знак (*), то данная модель может привести к плохой специфичности метода причинности по Грейнджеру. Если рядом с полученным числом стоит знак (o), то данная модель может привести к плохой чувствительности метода причинности по Грейнджеру. At -шаг дискретизации, равный 1/512 секунд.
**If a sign (*) is placed next to the obtained number, this model may lead to poor specificity of Granger causality method. If a sign (o) is placed next to the obtained number, this model may lead to poor sensitivity of the Granger causality method. At - sampling step equal to 1/512 s.
ложить, что это отношение должно быть на пару порядков больше, тогда результаты, полученные на системах ФитцХью-Нагумо и ван дер Поля-Тоды, будут очень похожи на результаты по экспериментальным данным.
2.2. Анализ связанности для тестовых систем. Для каждого «разряда» тестовых систем при исследовании причинности по Грейнджеру были использованы 30-секундные интервалы: 10 с до начала «разряда», 10 с «разряда» и 10 с после окончания «разряда». На графиках зависимости улучшения прогноза Р1теап от времени Ь чёрные сплошные вертикальные линии показывают начало и конец «разряда», расстояние от чёрных штриховых вертикальных линий до чёрных сплошных линий показывает ширину скользящего окна. Если изменения улучшения прогноза начались в пределах окна (от штриховой до сплошной линии), то они могут быть обусловлены тем, что окно захватывает переходной процесс от «фона» к «разряду». Соответственно, изменения улучшения прогноза уже не относятся к преиктальной активности. Здесь может наблюдаться феномен «уши». Красные точки соответствуют значению Р1 значимо большему фонового уровня на уровне значимости 0.05, а синие
точки - значимо меньшему. Над каждым графиком написано влияние какой системы на какую проверяется у ^ х.
Анализ связанности методом причинности по Грейнджер для тестовых систем показал следующие результаты. Для всех систем при единичной дальности прогноза метод выявляет ложные связи (плохая специфичность). Подобную картину можно видеть на примере системы Хиндмарш-Розе на рис. 3. Причём плохая специфичность оказывается при анализе в скользящем окне шириной и и> = 1 с, и и> = 2 с. Для системы Хиндмарш-Розе при малых дальностях прогноза метод детектирует двунаправленное взаимодействие между первым и вторым осциллятором вместо однонаправленного воздействия со стороны первого осциллятора на второй. А для систем ФитцХью-Нагумо и Хиндмарш-Розе ложные связи выявляются ещё и со стороны первого и второго осциллятора на взаимосвязанные третий и четвёр-
1.0 0.8 0.6 0.4 0.2 0.0,
1 1->2
: :
• ;
1.0 0.8 0.6 0.4 0.2 0.0,
1->3
:
*
I* ; 1
! :
1.0 0.8 0.6 0.4 0.2 0.0,
1-И
1.0 0.8 0.6 0.4 0.2 0.0,
4-1
J
0 5 10 15 20 25 30
0 5 10 15 20 25 30
3^-1
|
• :
0 5 10 15 20 25 30 0 5 10 15 20 25 30 1.0 0.8 0.6 0.4 0.2 0.0,
I 3->2
:
; 1 J
• W ;.
г Jf
0 5 10 15 20 25 30
\ ?->3
ь
р 0 «
0 5 10 15 20 25 3
2^4
0 5 10 1 5 20 25 30 0 5
¿4 '
10 15 20 25 30
1.0 0.8 0.6 0.4 0.2 0.0,
: 3-И
■ 1 ' A j ■
| 4^1
;
!» t
I.
А /
0 5 1 0 15 2 0 25 3
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
4->2|
; ■
%
/
0 5 10 15 20 25 30
4->3
•
*
0 5 10 15 20 25 30
Рис. 3. Результат статистического анализа полученных значений улучшения прогноза для связанных систем Хиндмарш-Розе при т = 1, w = 2 с. По оси абсцисс отложено время в секундах. Начало и конец «разряда» отмечены чёрными вертикальными линиями. Расстояние между чёрной вертикальной штриховой и сплошной линией показывает ширину скользящего окна. По оси ординат отложено улучшение прогноза PImean, чёрные точки соответствуют значению PImean значимо большему значения фонового уровня PIbg на уровне значимости 0.05. Над каждым графиком написано влияние какой системы на какую проверяется y ^ x
Fig. 3. The statistical analysis of the obtained prediction improvement values for the coupled Hindmarsh-Rose systems at т = 1, w = 2 s. X-axis: time in seconds. X-axis: time in seconds. The «seizure» onset and the «seizure» termination are marked with the black vertical lines. The distance between the black vertical dashed line and the solid line shows the moving window length. Y-axis: the prediction improvement PImean, black points indicate values PImean significantly larger than baseline level PIbg at the significance level of 0.05. In the title of each graph driving and driven systems are marked y ^ x
тый осцилляторы. Полученный результат хорошо согласуется с данными таблицы 4. Это распространённая ситуация, описанная в [29], возникает в случае синхронных колебаний двух систем. То есть модель с единичной дальностью прогноза помогает определять не реальное взаимодействие систем, а всего лишь синхронность их колебаний. Метод детектирует ложноположительные связи не только при единичной дальности прогноза, но и при близких к ней. При больших дальностях прогноза, начиная с т = Т/32, специфичность метода становится хорошей. Похожий результат был получен в [23] для систем ФитцХью-Нагумо, Рёсслера и Лоренца, а в данной работе дополнен результатами на системах Хиндмарш-Розе и ван дер Поля-Тоды.
При дальности прогноза т = Т/16 при w = 1 с специфичность не очень хорошая (рис. 4), но при w = 2 си специфичность, и чувствительность оптимальные. Стоит отметить, что для больших дальностей прогноза т = Т/12, Т/8 и Т/4 графики
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
1—>2
.......... ........ щ ......... Ä......
I
; I44!
;
t....... ......... ......
I I
1.0 0.8 0.6 0.4 0.2 0.0,
2-М
0 5 10 15 20 25 30
0 5 10 15 20 25 30
0 5 10 15 20 25 30
+3:
'.........I.......! m '1 й.......1...........
0 5 0 1 5 J 0 25 3
0 5 10 15 20 25 30
1.0 0.8 0.6 0.4 0.2 0.0,
?- 44
I ijj ..........
: a*
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
+1 i
.......... ........ .......... f Иг
i i to ..........
;
0 5 10 15 20 25 30
1.0 0.8 0.6 0.4 0.2 0.0,
344
...........I........! .......A.....*
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
4- ■+1 i
......... ........ .....' ? .......... ..........
; [ 4->2
..........I"......! I *
44-3!
.....a if" .......'
(JN
;
0 5 10 15 20 25 30
Рис. 4. Результат статистического анализа полученных значений улучшения прогноза для связанных систем ФитцХью-Нагумо при т = Т/16, w = 1 с. По оси абсцисс отложено время в секундах. Начало и конец «разряда» отмечены чёрными вертикальными линиями. Расстояние между чёрной вертикальной штриховой и сплошной линией показывает ширину скользящего окна. По оси ординат отложено улучшение прогноза PImean, чёрные точки соответствуют значению PImean значимо большему значения фонового уровня PIbg на уровне значимости 0.05. Над каждым графиком написано влияние какой системы на какую проверяется y ^ x
Fig. 4. The statistical analysis of the obtained prediction improvement values for the coupled FitzHugh-Nagumo systems at т = t/16, w = 1 s. X-axis: time in seconds. X-axis: time in seconds. The «seizure» onset and the «seizure» termination are marked with the black vertical lines. The distance between the black vertical dashed line and the solid line shows the moving window length. Y-axis: the prediction improvement PImean, black points indicate values PImean significantly larger than baseline level PIbg at the significance level of 0.05. In the title of each graph driving and driven systems are marked y ^ x
качественно не отличаются, и для -ш = 1 с все связи определяются правильно, только график выглядит менее «гладким», чем для и> = 2 с. Из чего можно сделать вывод, что чем меньше дальность прогноза, тем сильнее результат зависит от эффективной длины ряда, по которой строится математическая прогностическая модель. Так как для нас крайне важно как можно лучшее временное разрешение полученных результатов, то есть смысл взять как можно большую из возможных дальностей прогноза, чтобы использовать наименьшую из возможных ширину скользящего окна.
С целью проверки гипотезы о том, что даже неоптимальная с точки зрения критерия Шварца модель может подойти для метода причинности по Грейн-джеру [24], была построена линейная модель (Р = 1) при дальности прогноза т=Т/12. Для системы ФитцХью-Нагумо (рис. 5) результат оказался очень хорошим: детектированы все имеющиеся связи, лишних связей не обнаружено. Для системы
1-V?
к 0 1 нй 5
0 i «Л 1 i ? XV» 0 ? и» 5 3
0.4
1.0 0.8 0.6 0.4 0.2 0.0,
0 5 10 15 20 25 3
1->4
1.0 0.8 0.6 0.4 0.2 0.0,
24 i
' i1 :
1.0 0.8 0.6 0.4 0.2 0.0,
3
:
! i
3- +1
МАЙ \J>
0 1 0 15 2 0 25 3
1.0 0.8 0.6 0.4 0.2 0.0,
0 5 10 15 20 25 30
0 5 10 1 5 20 25 30 0 5 10 1 5 20 25 30
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
4-И
0 5 10 15 20 25 30
4-v?
4—
0 5 10 15 20 25 30 0 5 10 15 20 25 30 1.0 0.8 0.6 0.4 0.2 0.0,
.......игЩ
0 5 10 15 20 25 30
f
: I
Рис. 5. Результат статистического анализа полученных значений улучшения прогноза для связанных систем ван дер Поля-Тоды т = T/12, w = 2 с. По оси абсцисс отложено время в секундах. Начало и конец «разряда» отмечены чёрными вертикальными линиями. Расстояние между чёрной вертикальной штриховой и сплошной линией показывает ширину скользящего окна. По оси ординат отложено улучшение прогноза PImean, чёрные точки соответствуют значению PImean значимо большему значения фонового уровня PIbg на уровне значимости 0.05; синие точки - значимо меньшему. Над каждым графиком написано влияние какой системы на какую проверяется y ^ x
Fig. 3. The statistical analysis of the obtained prediction improvement values for for the coupled van der Pol-Toda systems т = T/12, w = 2 s. X-axis: time in seconds. X-axis: time in seconds. The «seizure» onset and the «seizure» termination are marked with the black vertical lines. The distance between the black vertical dashed line and the solid line shows the moving window length. Y-axis: the prediction improvement PImean, black points indicate values PImean significantly larger than baseline level PIbg at the significance level of 0.05. In the title of each graph driving and driven systems are marked y ^ x
Хиндмарш-Розе использование линейной модели приводит к потери чувствительности, что нетипично, так как обычно чем меньше коэффициентов у модели, тем лучше чувствительность метода, правда, и специфичность при этом ухудшается. Для системы ван дер Поля-Тоды линейная модель неплохо работает для случая двустороннего взаимодействия. В случае односторонней связи метод демонстрирует плохую специфичность (при исследовании воздействия с ведомой на ведущую систему 2 ^ 1 во время «разряда» появляются черные точки) и не очень хорошую чувствительность (при исследовании воздействия с ведущей на ведомую систему 1 ^ 2 во время «разряда» черные точки далеко не все), чем при использовании кубического полинома. Из всего вышесказанного о линейной модели можно сделать следующий вывод: для хорошей специфичности метода гораздо важнее выбрать правильную дальность прогноза, нежели подобрать правильную степень полинома.
По результатам анализа наблюдений для всех трёх систем следует отметить некоторые специфические тонкости параметризации:
1. При единичном лаге часто получается плохо обусловленная матрица, поэтому метод причинности по Грейнджеру начинает выдавать значения, выходящие за диапазон [0; 1].
2. При порядке полинома Р = 6 (для скользящего окна в 2 с именно такой порядок по критерию Шварцу подбирается без дополнительного ограничения на совместную модель) получается высокое значение Р1 даже в фоне (около 0.5), падает чувствительность (не видит связи там, где она есть).
3. Аналогично чувствительность падает для Р=4 и ширине скользящего окна 1 с. Перечисленные тонкости процесса параметризации подтверждают правильность выбранного нами более сильного ограничения на выбор эффективной длины ряда для моделирования.
2.3. Анализ связанности по сигналам экспериментальных внутричерепных ЭЭГ. Для каждого разряда экспериментальных внутричерепных ЭЭГ при исследовании причинности по Грейнджеру были использованы 40-секундные интервалы: 10 с до начала разряда, 10 с от начала разряда, 10 с перед окончанием разряда и 10 с после окончания разряда. Полученные зависимости улучшения прогноза от времени Р1 (¿) усреднялись по всем разрядам путём совмещения начал или концов, и получалась зависимость Р1теап (¿), которая и откладывалась на графиках.
Для экспериментальных данных ЭЭГ при использовании единичной дальности прогноза во время разряда детектируется усиление влияния со стороны затылочной коры на гиппокамп, но затем наблюдается чёткая двусторонняя связь (рис. 6, а), начинающаяся в тот момент, когда на спектрограмме пропадают высокие частоты и падает амплитуда (см. рис. 1). Скорее всего, в этот момент реальное взаимодействие уже прекратилось, но колебания, генерируемые обеими структурами мозга, всё ещё остаются синхронными. Но как было показано на тестовых примерах, метод причинности по Грейнджеру при малых дальностях прогноза демонстрирует плохую специфичность. Следовательно, не стоит доверять результатам, полученным с использованием подобной модели, и на экспериментальных данных.
В случае использования линейной модели при дальности прогноза т = Т/12 также можно видеть продолжающееся увеличенное по сравнению с фоном воздействие со стороны коры на гиппокамп после окончания разряда (рис. 6, Ь). Как видно
на системе ван дер Поля-Тоды, линейная модель действительно обладает не очень хорошей специфичностью и может детектировать несуществующую связь.
Для других дальностей прогноза т = Т/16, T/12, Т/8, Т/4 при использовании квадратичного полинома двусторонняя связь не наблюдается нигде (рис. 6, c, d). Видно только усиление воздействия со стороны затылочной коры на гиппокамп почти в пределах разряда (между чёрными вертикальными линиями): значимое усиление воздействия со стороны затылочной коры на гиппокамп начинается примерно за 2 с до размеченного совместно с физиологами начала разряда и заканчивается ровно в момент размеченного конца разряда. При ширине скользящего окна w = 1 с получаются такие же выводы об увеличении связанности, что и при w = 2 с. Исходя из результатов, полученных на модельных системах, в этом случае можно говорить о возможности использования такого узкого скользящего окна (4 характерных пе-
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 о.о,
10
10
Нр-ЮС
15
20
ОС->Нр
15
20
a
25
30
35
40
25 30 35 40
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
10
10
Нр-ЮС
15
20
25
30
35
40
ОС-»Нр
/
•i-AA'
15 20 25 30 35 40
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
.....Hp-ioc
.................. 1 14 й i iW*^
\ ОС-»Нр :
;
10
15
20
С
25
30
35
40
1.0 0.8 0.6 0.4 0.2 0.0,
1.0 0.8 0.6 0.4 0.2 0.0,
10
Hp- vnc,
SM.
ОС->Нр
15
20
d
25
30
35
40
Рис. 6. Результат статистического анализа полученных значений улучшения прогноза для реальных лимбических разрядов: а - т = 1, w = 2 с, P = 2; b - т = T/12,w = 2 с, P = 1; c - т = T/12, w = 1 с, P = 2; d - т = T/12, w = 2 с, P = 2. По оси абсцисс отложено время в секундах. Начало и конец разряда отмечены чёрными вертикальными линиями. Расстояние между чёрной вертикальной штриховой и сплошной линией показывает ширину скользящего окна. По оси ординат отложено улучшение прогноза PImean, чёрные точки соответствуют значению PImean значимо большему фонового уровня PIbg на уровне значимости 0.05
Fig. 6. The statistical analysis of the obtained prediction improvement values (a) for the coupled FitzHugh-Nagumo systems т = t/16, w = 1 s, (b) for the coupled Hindmarsh-Rose systems т = 1, w = 2 s, (c) for the coupled van der Pol-Toda systems т = T/12, w = 2 s. X-axis: time in seconds. X-axis: time in seconds. The «seizure» onset and the «seizure» termination are marked with the black vertical lines. The distance between the black vertical dashed line and the solid line shows the moving window length. Y-axis: the prediction improvement PImean, black points indicate values PImean significantly larger than baseline level PIbg at the significance level of 0.05. In the title of each graph driving and driven systems are marked y ^ x
5
5
5
5
b
5
5
риода) для лучшего временного разрешения. Но если бы результаты, полученные с использованием более широкого окна, оказались отличными от результатов для узкого окна, то такую возможность нельзя было использовать.
Заключение
В настоящее время накопилось значительное число результатов применения методов оценки связанности (включая метод причинности по Грейнджеру) к задаче детектирования взаимодействия структур мозга, в том числе при эпилепсии (см., например, [5, 12, 13,15,21]). При этом часто надёжность выводов оказывается под вопросом из-за несовершенства используемых моделей, чему посвящён ряд специальных работ, в частности, [29, 35]. Верная интерпретация результатов обычно требует значительного времени и дополнительных исследований, как случилось с выводами работы [21], переосмысленными в [25]. Основным путём к повышению чувствительности и специфичности в настоящее время является тщательный подбор структуры модели и вида нелинейных функций, как это было сделано в [11,38]. Однако в [21] подбор размерностей моделей был осуществлён, но это не помешало получить неинформативные результаты, затем неверно истолкованные. Следовательно, использование только универсальных критериев вроде критерия Шварца [39] для подбора параметров явно недостаточно.
Идея настоящей работы заключается в том, что повысить надёжность выводов о связанности можно, апробируя методы на симулированных данных, полученных от моделей исследуемых сигналов, воспроизводящих значительное число свойств экспериментальных временных рядов, например, таких как модели, предложенные для абсансной эпилепсии в работах [42,43]. Хотя такой подход более трудоёмкий и до некоторой степени граничит с искусством [1], он обещает лучшее понимание и более надёжную интерпретацию полученных значений мер связанности, что часто важнее увеличения числа исследованных пациентов, животных или здоровых испытуемых, или использования более мощных методов статистического тестирования на значимость.
Благодаря использованию модельных систем для тестирования метода причинности по Грейнджеру, можно сделать следующий довольно надёжный вывод: во время вызванных введением агониста эндоканнабиноидных рецепторов лимбиче-ских разрядов примерно за 2 с до начала приступа наблюдается значимое усиление воздействия со стороны затылочной коры на гиппокамп. Данное воздействие ослабевает до фонового уровня ровно в момент конца разряда. Важность данного вывода подчёркивается тем, сколько усилий в настоящее время прилагается для изучения взаимодействий и локализации эпилептического фокуса при вторично генерализованных формах эпилепсии [44].
Авторы выражают особую благодарность Илье Вячеславовичу Сысоеву, доценту Саратовского государственного университета, за обсуждения результатов и ценные советы, а также Tineke van Rijn, assistant professor at Radboud University, Nijmegen, the Netherlands, за предоставленные экспериментальные данные. Работа выполнена при поддержке Стипендии Президента для поддержки молодых учёных СП-3605.2018.4 и РФФИ, грант 17-02-00307 и грант 18-015-00418.
Библиографический список
1. Bezruchko B.P., Smirnov D.A. Extracting Knowledge From Time Series. Berlin: Springer, 2010.
2. Granger C.W.J. Investigating causal relations by econometric models and cross-spectral methods // Econometrica. 1969. Vol. 37(3). P. 424-438.
3. van Rijn C., Gaetani S., Santolini I., Badura A., Gabova A., Fu J., Watanabe M., Cuomo V., van Luijtelaar G., Nicoletti F., Ngomba R. WAG/Rij rats show a reduced expression of CB1 receptors in thalamic nuclei and respond to the CB1 receptor agonist, R(+)WIN55,212-2, with a reduced incidence of spike-wave discharges // Epilepsia. 2010. Vol. 51(8). P. 1511-1521.
4. Liittjohann A., van Luijtelaar G. The dynamics of cortico-thalamo-cortical interactions at the transition from pre-ictal to ictal LFPs in absence epilepsy // Neurobiology of Disease. 2012. Vol. 47. P. 47-60.
5. Liittjohann A., Schoffelen J., van Luijtelaar G. Termination of ongoing spike-wave discharges investigated by cortico-thalamic network analyses // Neurobiology of Disease. 2014. Vol. 70. P. 127-137.
6. Колосов А.В., Нуйдель И.В., Яхно В.Г. Исследование динамических режимов в математической модели элементарной таламокортикальной ячейки // Известия вузов. Прикладная нелинейная динамика. 2016. Т. 24(5). С. 72-83.
7. Blumenfeld H., Varghese G., Purcaro M., Motelow J., Enev M., McNally K., Levin A., Hirsch L., Tikofsky R., Zubal I., Paige A., Spencer S. Cortical and subcortical networks in human secondarily generalized tonic-clonic seizures // Brain. 2009. Vol. 1324. P. 999-1012.
8. Wallace M., Blair R., Falenski K., Martin B., De Lorenzo R. The endogeneous cannabinoid system regulates seizure frequency and duration in a model of temporal lobe epilepsy // J Pharmacol. Exp. Ther. 2003. Vol. 307. P. 129-137.
9. HaneefZ., Lenartowicz A., Yeh H., Levin H., Engel J., Stern J.Functional connectivity of hippocampal networks in temporal lobe epilepsy // Epilepsia. 2014. Vol. 551. P. 137-145.
10. Ding M., Chen Y., Bressler S. Granger causality: basic theory and application to neuroscience - Handbook of time series analysis. In: «Handbook of Time Series Analysis: Recent Theoretical Developments and Applications», edited by Bjorn Schelter, Matthias Winterhalder, Jens Timmer, 2006, Wiley-VCH Verlag GmbH & Co. KGaA.
11. Chen Y., Rangarajan G., Feng J., Ding M. Analyzing Multiple Nonlinear Time Series with Extended Granger Causality // Physics Letters A. 2004. Vol. 324(1). P. 26-35.
12. Marinazzo D., Pellicoro M., Stramaglia S. Nonlinear parametric model for Granger causality of time series // Phys. Rev. E. 2006. Vol. 73. 066216.
13. Marinazzo D., Pellicoro M., Stramaglia S. Kernel Method for Nonlinear Granger Causality // Phys. Rev. Lett. 2008. Vol. 100. 144103.
14. Lehnertz K., Andrzejak R., Arnhold J., Kreuz T., Mormann F., Rieke C., Widman G., Elger C. Nonlinear EEG Analysis in Epilepsy: Its Possible Use for Interictal Focus Localization, Seizure Anticipation and Prevention // Journal of Clinical Neurophysiology. 2001. Vol. 18(3). P. 209-222.
15. Gourevitch B., Le Bouquin-Jeanrns R., Faucon G. Linear and nonlinear causality
between signals: methods, examples and neurophysiological applications // Biological Cybernetics. 2006. Vol. 95. P. 349-369.
16. Pereda E., QuianQuiroga R., Bhattacharya J. Nonlinear multivariate analysis of neurophysiological signals // Progress in Neurobiology. 2005. Vol. 77(1). P. 1-37.
17. Cekic S., Grandjean D., Renaud O. Time, frequency and time-varying causality measures inNeuroscience// Statistics inMedicine. 2018. Vol. 37(11). P. 1910-1931.
18. Sysoeva M., Luttjohann A., van Luijtelaar G., Sysoev I. Dynamics of directional coupling underlying spike-wave discharges//Neuroscience. 2016. Vol. 314. P.75-89.
19. Сысоева М.В., Ситникова Е.Ю., Сысоев И.В. Таламо-кортикальные механизмы инициации, поддержания и прекращения пик-волновых разрядов у крыс WAG/Rij // Журнал высшей нервной деятельности имени И.П. Павлова. 2016. Т. 66(1). С. 103-112.
20. Sysoeva M.V., Vinogradova L.V., Kuznetsova G.D., Sysoev I.V., van Rijn C.M. Changes in cortico-cortical and cortico-hippocampal network during absence seizures in WAG/Rij rats revealed with time varying Granger causality // Epilepsy & Behavior. 2016. Vol. 64. P. 44-50.
21. Sitnikova E., Dikanev T., Smirnov D., Bezruchko B., van Luijtelaar G. Granger causality: Cortico-thalamic interdependencies during absence seizures in WAG/Rij rats // Journal of Neuroscience Methods. 2008. Vol. 170. P. 245-254.
22. Sysoeva M.V., Sysoev I.V. Mathematical modeling of encephalogram dynamics during epileptic seizure // Technical Physics Letters. 2012. Vol. 38(2). P. 151-154.
23. Сысоева М.В., Диканев Т.В., Сысоев И.В. Выбор временных масштабов при построении эмпирической модели // Изв. вузов. Прикладная нелинейная динамика. 2012. Т. 20(2). С. 54-62.
24. Корнилов М.В., Сысоев И.В. Влияние выбора структуры модели на работоспособность метода нелинейной причинности по Грейнджеру // Изв. вузов. Прикладная нелинейная динамика. 2013. T. 21(2). С. 3-16.
25. Sysoeva M.V., Sitnikova E., Sysoev I.V., Bezruchko B.P., van Luijtelaar G. Application of adaptive nonlinear Granger causality: Disclosing network changes before and after absence seizure onset in a genetic rat model // Journal of Neuroscience Methods. 2014. Vol. 226. P. 33-41.
26. Zou C., Feng J.Granger causality vs. dynamic Bayesian network inference: a comparative study // BMC Bioinformatics. 2009. Vol. 10: 122.
27. Prokhorov M.D., Ponomarenko V.I. Estimation of coupling between time-delay systems from time series // Phys. Rev. E. 2005. Vol. 72. 016210.
28. Sysoev I.V., Prokhorov M.D., Ponomarenko V.I., Bezruchko B.P. Reconstruction of ensembles of coupled time-delay system from time series // Phys. Rev. E. 2014. Vol. 89. 062911.
29. Sysoev I.V., Sysoeva M.V. Detecting changes in coupling with Granger causality method from time series with fast transient processes // Physica D: Nonlinear Phenomena. 2015. Vol. 309. P. 9-19.
30. Безручко Б.П., Смирнов Д.А., Зборовский А.В., Сидак Е.В., Иванов Р.Н., Беспя-тов А.Б. Реконструкция по временному ряду и задачи диагностики // Технологии живых систем. 2007. Т. 4(3). С. 49-56.
31. Besruchko B.P., Smirnov D.A. Constructing nonautonomous differential equations from experimental time series // Phys. Rev. E. 2000. Vol. 63. 016207.
32. Безручко Б.П., Смирнов Д.А., Сысоев И.В., Селезнев Е.П. Реконструкция моделей неавтономных систем с дискретным спектром воздействия // Письма в ЖТФ. 2003. Т. 29(19), С. 69-76.
33. Packard N., Crutchfield J., Farmer J., Shaw R. Geometry from a Time Series // Phys. Rev. Lett. 1980. Vol. 45. P. 712-716.
34. Legendre A. Appendice sur la methodes des moindres quarres. Nouvelles methodes pour la determination des orbites des cometes // Paris: Firmin-Didot. 1805. (in French)
35. Smirnov D.A, Bezruchko B.P. Spurious causalities due to low temporal resolution: Towards detection of bidirectional coupling from time series // Europhys. Lett. 2012. Vol. 100. 10005.
36. Smirnov D.A. Mokhov I.I. From Granger causality to long-term causality: Application to climatic data // Physical Review E. 2009. Vol.80. 016208.
37. Hesse W., Moller E., Arnold M., Schack B. The use of time-variant EEG Granger causality for inspecting directed interdependencies of neural assemblies // Journal of Neuroscience Methods. 2003. Vol. 124. P. 27-44.
38. Kornilov M.V., Medvedeva T.M., Bezruchko B.P., Sysoev I.V. Choosing the optimal model parameters for Granger causality in application to time series with main timescale // Chaos, Solitons & Fractals. 2016. Vol. 82. P. 11-21.
39. Schwarz G. Estimating the Dimension of a Model // The Annals of Statistics. 1978. Vol. 6(2). P. 461-464.
40. Rosenstein M., Collins J., De Luca C. A practical method for calculating largest Lyapunov exponents from small data sets // Physica D. 1993. Vol. 6. P. 117-134.
41. Student. The probable error of a mean. Biometrika. 1908. 6(1). P. 1-25.
42. Suffczynski P., Kalitzin S., Lopes da Silva F.H. Dynamic sofnon-convulsive epileptic phenomena modeled by a bistable neuronal network // Neuroscience. 2004. Vol. 126. P. 467-484.
43. Сычева М.В., Кузнецова Г.Д., Сы^ев И.В. Модел^ование cигналов элекфо-энцефало^амм ^ью ^и абcанcной эпилепоти в ^иложении к анализу cвя-занноети между отделами мозга// Биофизика. 2016. Т. 61. В. 4. С. 782-792.
44. Jeanne T. Paz J.T., Huguenard J.R. Microcircuits and their interactions in epilepsy: is the focus out of focus? //Nature Neuroscience. 2015. Vol. 18, iss. 3. Pp. 351-359. Doi: 10.1038/nn.3950.
References
1. Bezruchko B.P., Smirnov D.A. Extracting Knowledge From Time Series. Berlin: Springer, 2010.
2. Granger C.W.J. Investigating Causal Relations by Econometric Models and Cross-Spectral Methods. Econometrica, 1969, vol. 37(3), p. 424-438.
3. van Rijn C., Gaetani S., Santolini I., Badura A., Gabova A., Fu J., Watanabe M., Cuomo V., van Luijtelaar G., Nicoletti F., Ngomba R. WAG/Rij rats show a reduced expression of CB1 receptors in thalamic nuclei and respond to the CB1 receptor agonist, R(+)WIN55,212-2, with a reduced incidence of spike-wave discharges. Epilepsia, 2010, vol. 51(8), p. 1511-1521.
4. Luttjohann A., van Luijtelaar G. The dynamics of cortico-thalamo-cortical interactions at the transition from pre-ictal to ictal LFPs in absence epilepsy. Neurobiology of Disease, 2012, vol. 47, p. 47-60.
5. Luttjohann A., Schoffelen J., van Luijtelaar G. Termination of ongoing spike-wave discharges investigated by cortico-thalamic network analyses // Neurobiology of Disease, 2014, vol. 70, p. 127-137.
6. Kolosov A.V., Nuidel I.V., Yakhno V.G. Research of dynamic modes in the mathematical model of elementary thalamocortical cell // Izvestiya VUZ. Applied Nonlinear Dynamics, 2016, vol. 24(5), p. 72-83. (in Russian)
7. Blumenfeld H., Varghese G., Purcaro M., Motelow J., Enev M., McNally K., Levin A., Hirsch L., Tikofsky R., Zubal I., Paige A., Spencer S. Cortical and subcortical networks in human secondarily generalized tonic-clonic seizures. Brain, 2009, vol. 1324, p. 999-1012.
8. Wallace M., Blair R., Falenski K., Martin B., De Lorenzo R. The endogeneous cannabinoid system regulates seizure frequency and duration in a model of temporal lobe epilepsy. J Pharmacol. Exp. Ther., 2003, vol. 307, p. 129-137.
9. Haneef Z., Lenartowicz A., Yeh H., Levin H., Engel J., Stern J. Functional connectivity of hippocampal networks in temporal lobe epilepsy. Epilepsia, 2014, vol. 551, p. 137-145.
10. Ding M., Chen Y., Bressler S. Granger causality: basic theory and application to neuroscience - Handbook of time series analysis. In: «Handbook of Time Series Analysis: Recent Theoretical Developments and Applications», edited by Bjorn Schelter, Matthias Winterhalder, Jens Timmer, 2006, Wiley-VCH Verlag GmbH & Co. KGaA.
11. Chen Y., Rangarajan G., Feng J., Ding M. Analyzing Multiple Nonlinear Time Series with Extended Granger Causality. Physics Letters A., 2004, vol. 324(1), p. 26-35.
12. Marinazzo D., Pellicoro M., Stramaglia S. Nonlinear parametric model for Granger causality of time series. Phys. Rev. E., 2006, vol. 73, 066216.
13. Marinazzo D., Pellicoro M., Stramaglia S. Kernel Method for Nonlinear Granger Causality. Phys. Rev. Lett., 2008, vol. 100, 144103.
14. Lehnertz K., Andrzejak R., Arnhold J., Kreuz T., Mormann F., Rieke C., Widman G., Elger C. Nonlinear EEG Analysis in Epilepsy: Its Possible Use for Interictal Focus Localization, Seizure Anticipation and Prevention. Journal of Clinical Neurophysiology, 2001, vol. 18(3), p. 209-222.
15. Gourevitch B., Le Bouquin-Jeannes R., Faucon G. Linear and nonlinear causality between signals: methods, examples and neurophysiological applications. Biological Cybernetics, 2006, vol. 95, p. 349-369.
16. Pereda E., QuianQuiroga R., Bhattacharya J. Nonlinear multivariate analysis of neurophysiological signals. Progress in Neurobiology, 2005, vol. 77(1), p. 1-37.
17. Cekic S., Grandjean D., Renaud O. Time, frequency and time-varying causality measures in Neuroscience. Statistics in Medicine, 2018, vol.37(11), p. 1910-1931.
18. Sysoeva M., Luttjohann A., van Luijtelaar G., Sysoev I. Dynamics of directional coupling underlying spike-wave discharges. Neuroscience, 2016, vol. 314, p. 75-89.
19. Sysoeva M.V., Sitnikova E., Sysoev I.V. Thalamo-cortical mechanisms of initiation, maintenance and termination of spike-wave discharges at WAG/Rij rats. Zhurnal vysshei nervnoi deyatel'nosti im. I.P. Pavlova, 2016, vol. 66(1), p. 103-112. (in Russian)
20. Sysoeva M.V., Vinogradova L.V., Kuznetsova G.D., Sysoev I.V., van Rijn C.M. Changes in cortico-cortical and cortico-hippocampal network during absence seizures
in WAG/Rij rats revealed with time varying Granger causality. Epilepsy & Behavior, 2016, vol. 64, p. 44-50.
21. Sitnikova E., Dikanev T., Smirnov D., Bezruchko B., van Luijtelaar G. Granger causality: Cortico-thalamic interdependencies during absence seizures in WAG/Rij rats. Journal of Neuroscience Methods, 2008, vol. 170, p. 245-254.
22. Sysoeva M.V., Sysoev I.V. Mathematical modeling of encephalogram dynamics during epileptic seizure. Technical Physics Letters, 2012, vol. 38(2), p. 151-154.
23. Sysoeva M.V., Dikanev T.V. and Sysoev I.V. Selecting time scales for empirical model construction. Izvestiya VUZ, Applied Nonlinear Dynamics, 2012, vol. 20(2), p. 54-62. (in Russian)
24. Kornilov M.V., Sysoev I.V. Influence of the choice of the model structure for working capacity of nonlinear Granger causality approach. Izvestiya VUZ, Applied Nonlinear Dynamics, 2013, vol. 21(2), p. 3-16. (in Russian)
25. Sysoeva M.V., Sitnikova E., Sysoev I.V., Bezruchko B.P., van Luijtelaar G. Application of adaptive nonlinear Granger causality: Disclosing network changes before and after absence seizure onset in a genetic rat model. Journal ofNeuroscience Methods, 2014, vol. 226, p. 33-41.
26. Zou C., Feng J. Granger causality vs. dynamic Bayesian network inference: a comparative study. BMC Bioinformatics, 2009, vol. 10: 122.
27. Prokhorov M.D., Ponomarenko V.I. Estimation of coupling between time-delay systems from time series. Phys. Rev. E., 2005, vol. 72, 016210.
28. Sysoev I.V., Prokhorov M.D., Ponomarenko V.I., Bezruchko B.P. Reconstruction of ensembles of coupled time-delay system from time series. Phys. Rev. E., 2014, vol. 89, 062911.
29. Sysoev I.V., Sysoeva M.V. Detecting changes in coupling with Granger causality method from time series with fast transient processes. Physica D: Nonlinear Phenomena, 2015, vol. 309, p. 9-19.
30. Bezruchko B.P., Smirnov D.A., Zborovsky A.B., Sidak E.V., Ivanov R.N., Bespyatov A.B. Reconstruction of the time series and diagnostic problems. Technologies of living systems, 2007, vol. 4(3), p. 49-56. (in Russian)
31. Besruchko B.P., Smirnov D.A. Constructing nonautonomous differential equations from experimental time series. Phys. Rev. E., 2000, vol. 63, 016207.
32. Smirnov D.A., Sysoev I.V., Seleznev Ye.P., Bezruchko B.P. Reconstructing nonauto-nomous system models with discrete spectrum of external action. Technical Physics Letters, 2003, vol. 29(10), p. 824-827.
33. Packard N., Crutchfield J., Farmer J., Shaw R. Geometry from a Time Series. Phys. Rev. Lett., 1980, vol. 45, p. 712-716.
34. Legendre A. Appendice sur la methodes des moindres quarres. Nouvelles methodes pour la determination des orbites des cometes. Paris: Firmin-Didot. 1805. (in French)
35. Smirnov D.A. Bezruchko B.P. Spurious causalities due to low temporal resolution: Towards detection of bidirectional coupling from time series. Europhys. Lett., 2012, vol. 100, 10005.
36. Smirnov D.A. Mokhov I.I. From Granger causality to long-term causality: Application to climatic data. Physical Review E., 2009, vol. 80, 016208.
37. Hesse W., Moller E., Arnold M., Schack B. The use of time-variant EEG Granger
causality for inspecting directed interdependencies of neural assemblies. Journal of Neuroscience Methods, 2003, vol. 124, p. 27-44.
38. Kornilov M.V., Medvedeva T.M., Bezruchko B.P., Sysoev I.V. Choosing the optimal model parameters for Granger causality in application to time series with main timescale. Chaos, Solitons & Fractals, 2016, vol. 82, p. 11-21.
39. Schwarz G. Estimating the Dimension of a Model. The Annals of Statistics. 1978, vol. 6(2), p. 461-464.
40. Rosenstein M., Collins J., De Luca C. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D., 1993, vol. 6, p. 117-134.
41. Student. The probable error of a mean. Biometrika, 1908, 6(1), p. 1-25.
42. Suffczynski P., Kalitzin S., Lopes da Silva F.H. Dynamic sofnon-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience, 2004, vol. 126, p. 467-484.
43. Sysoeva M.V., Kuznetsova G.D., Sysoev I.V. Modelling EEG signals from rats when analysing absence epilepsy in application to analysis of coupling between brain areas. Biophysics, 2016, vol. 61(4), p. 661-669.
44. Jeanne T. Paz, John R. Huguenard. Microcircuits and their interactions in epilepsy: is the focus out of focus? Nature Neuroscience, 2015, vol. 18, iss. 3, pp. 351-359. Doi: 10.1038/nn.3950.
Сысоева Марина Вячеславовна. Родилась в Саратове (1987). Окончила Лицей № 37 (2005, Саратов). Окончила с отличием бакалавриат (2009) и магистратуру (2011) факультета нано- и биомедицинских технологий СГУ имени Н.Г. Чернышевского по направлению «Биомедицинская инженерия». Научный руководитель - доцент, к.ф.-м.н., Т.В. Диканев. Защитила диссертацию на соискание учёной степени кандидата физико-математических наук на тему «Особенности реализации метода причинности по Грейнджеру для исследования электроэнцефалограмм при абсансной эпилепсии» по специальностям «Биофизика» и «Радиофизика» (2015, СГУ). Научные руководители: доцент, к.ф.-м.н. И.В. Сысоев и профессор, д.ф.-м.н. В.И. Пономаренко. С 2015 года работает на кафедре «Радиоэлектроника и телекоммуникации» Саратовского государственного технического университета имени Гагарина Ю.А. Научные интересы: анализ временных рядов, нейронаука, математическое моделирование.
Россия, 410054 Саратов, Политехническая, 77
Саратовский государственный технический университет имени Гагарина Ю.А. E-mail: bobrichek@mail.ru
Медведева Татьяна Михайловна - родилась в 1993 году в городе Энгельсе Саратовской области, в 2016 году окончила магистратуру по направлению «биотехнические системы и технологии» Саратовского государственного университета. С 2014 по 2017 год работала на кафедре динамического моделирования и биомедицинской инженерии СГУ, в 2016 году поступила в аспирантуру на эту кафедру. С 2018 года является младшим научным сотрудником Института высшей нервной деятельности и нейрофизиологии РАН. Автор 6 научных статей.
Россия, 117485 Москва, Бутлерова, 5А
Институт высшей нервной деятельности и нейрофизиологии РАН E-mail: Golova93tanya@gmail.com
