Роль нелинейности модели в диагностике связей при патологическом треморе методом грейнджеровской причинности Текст научной статьи по специальности «Математика»

Прикладные задачи

^^^^^^^^^^»нелинейной теории колебаний и вслн

Изв. вузов «ПНД», т. 18, № 4, 2010 УДК 530.182, 57.087.1

РОЛЬ НЕЛИНЕЙНОСТИ МОДЕЛИ В ДИАГНОСТИКЕ СВЯЗЕЙ ПРИ ПАТОЛОГИЧЕСКОМ ТРЕМОРЕ МЕТОДОМ ГРЕЙНДЖЕРОВСКОЙ ПРИЧИННОСТИ

И.В. Сысоев, А.С. Караваев, П.И. Наконечный

Оценка связи между системами различной природы - одно из наиболее востребованных направлений приложения методов нелинейной динамики. В данной работе сопоставляются возможности классического линейного метода оценки причинности по Грен-джеру и его нелинейного расширения на конкретных примерах как эталонных систем, так и на примере анализа нейрофизиологических данных. Показано, что нелинейный метод имеет большую чувствительность и чаще позволяет значимо детектировать связь.

Ключевые слова: Анализ связанности, нелинейная причинность по Гренджеру, обработка ЭЭГ, динамическое моделирование.

Введение

Развитие методов нелинейной динамики сопровождается их активным внедрением в область изучения объектов живой природы, что расширяет диагностический инструментарий медицины, физиологии, биологии. Возможности этих методик уже воспринимаются настроенными на новое врачами, как «добавление цвета к черно-белому изображению или звука к немому кино» [1, с. 93]. В частности, методы динамического моделирования востребованы для диагностики связей элементов сложных систем по дискретным записям их колебаний - временным рядам [2-4]. В статье рассматривается способ выявления таких связей с помощью относительного изменения точности прогноза поведения одного элемента при введении в прогностическую математическую модель данных о колебаниях другого элемента. При таком подходе уменьшение ошибки прогноза толкуется как признак воздействия второго элемента на первый. Эта мера, предложенная и реализованная с помощью авторегрессионных моделей К. Грейнджером [5], нобелевским лауреатом 2003 года по экономике, получила название «грейнджеровской причинности» и в настоящее время используется во многих областях знания.

Один из вариантов нелинейного расширения метода представлен в работе [6], где в качестве аппроксимирующих функций используются степенные полиномы общего вида. Возможность применения радиальных базисных функций для реализации нелинейной грейнджеровской причинности проиллюстрирована в [7], где такой подход применяется для анализа связанности частоты сердечных сокращений и давления крови в сосудах по экспериментальным данным. В работе [8] представлены результаты сравнения различных нелинейных расширений метода грейнджеровской причинности и других методов анализа, в частности, энтропии переноса. Также известен аналогичный метод анализа связанности сигналов по их фазам, предложенный в [9] и модернизированный в [10]; данная модификация применялась для анализа нейрофизиологических данных в [11]. Но, несмотря на ряд приведённых примеров, нелинейный вариант метода грейнджеровской причинности в настоящее время не получил широкого распространения, главным образом, в связи со сложностью выбора параметров метода.

Целью данной работы является сопоставление возможностей линейного и нелинейного методов грейнджеровской причинности на примере анализа нейрофизиологических данных с целью оценки связей между головным мозгом и конечностью при патологии cortical tremor. Выбор объекта обусловлен тем, что для этого вида патологии имеются физиологические основания априорно полагать наличие взаимной связи между областями моторной коры головного мозга и конечностью. Исходными для наших оценок являются временные ряды сигналов скальповой электроэнцефалограммы (ЭЭГ) и показаний акселерометра, закреплённого на руке. Достоверность результатов анализа реальной системы оценивалась путём контроля значимости с помощью суррогатных данных и реконструкции связей в системе связанных дифференциальных уравнений эталонных осцилляторов при различных соотношениях линейной и нелинейной компонент связи. Результаты анализа представленных примеров показывают, что при оценке грейнджеровской причинности упомянутые выше затраты на подбор параметров в нелинейной модели оправданы и неизбежны, поскольку нелинейный метод обладает большей чувствительностью в сравнении с линейным, позволяя детектировать связь в ситуациях, когда выявить её с помощью линейного подхода не удаётся.

В разделе 1 даётся подробное описание метода оценки грейнджеровской причинности, в разделе 2 с его помощью оцениваются связи между сигналами электрических потенциалов мозга и колебаний руки, а достоверность выводов и анализ результатов проведен в разделе 3 с помощью эталонной модели.

1. Описание метода

Основная идея метода была сформулирована в работе [5] для линейных моделей. Пусть у нас имеются два временных ряда: ряд {xn}N=1 от системы X и ряд [Vn}n=i от системы Y, где n - дискретное время, N - длина временного ряда. На основе анализа реализаций {xn}n=l и {уп}П=\ требуется определить, влияет ли система Y на систему X или нет.

На первом шаге строится индивидуальная модель

Xn = f (Xn-l,Xn-2l, Xn-Dsl, CS) + m (1)

где / - аппроксимирующая функция, I - лаг модели, - собственная размерность модели, с5 - неизвестные коэффициенты, а "п - остатки модели, по предположению, представляющие собою нормальный некоррелированный шум. В качестве функции / можно использовать разложение по некоторому базису, например, представить её как многомерный полином общего вида (см. [6]). Если имеется дополнительная информация об объекте X, то структуру функции / следует выбирать, основываясь на ней.

Коэффициенты с5 оцениваются методом наименьших квадратов по экспериментальной реализации {хп}1^=1. Полученная в результате модель имеет среднеквадратичную ошибку аппроксимации е\, равную дисперсии остатков

Следующим шагом строится совместная модель

Хп = д (хп-1 ,Хп-21, ■■■, Хп-Пе1,уп-1-%, Уп-21-т, ■■■, Уп-Оа 1-т, С) + , (2)

где для оценки коэффициентов аппроксимирующей функции д используется также реализация системы У; Оа размерность добавки (число учитываемых значений из ряда {уп}^=1); С - коэффициенты совместной модели; т - задержка, обусловленная конечным временем распространения сигнала (на практике часто неизвестна и должна быть определена).

Построив совместную модель (2) можно рассчитать среднеквадратичную ошибку прогноза ¿2. Случай е2 < е23 показывает, что данные из ряда системы У помогли предсказать поведение системы X. Тогда говорят, что У действует на X «по Грейн-джеру».

Следуя ряду работ, в частности [6], будем количественно характеризовать воздействие с помощью улучшения прогноза

Р1 =1 - ¿2. (3)

Р1 = 0 в случае, если данные из ряда У не помогают предсказывать динамику системы X, то есть. ¿2 = ¿1. Р1 достигает 1, если динамика X полностью описывается совместною моделью (е^ = 0), но не описывается индивидуальной.

Поскольку совместная ошибка может быть меньше индивидуальной вследствие действия случайных факторов, требуется оценка статистической значимости рассчитываемого значения Р1. Чтобы проверить, действительно ли введение зависимости от сигнала {уп}1^=1 в модель для сигнала {хп}1^=1 даёт качественное улучшение прогноза, не обусловленное случайными факторами, можно использовать различные подходы: получить оценку доверительного интервала из теоретических соображений о свойствах остатков модели, как в [9-11], либо использовать суррогатные временные ряды, как, например, в [12].

В данной работе используется анализ суррогатных данных, поскольку проверка статистической значимости не требует в этом случае предположений о свойствах исходных данных. Суррогатные данные готовятся из исходных сигналов путём случайного задания фаз компонент их фурье-образов [13].

Для оценки статистической значимости каждого рассчитываемого значения Р1 генерировался ансамбль из 100 рядов суррогатных данных и оценивался 95%-й квантиль. При расчётах Р1 (т), сопровождаемых перебором пробных значений задержки т, оценивался полный 95%-й квантиль РТаь3.

2. Анализ экспериментальных данных

2.1. Описание данных. Экспериментальная проверка работоспособности сопоставляемых методов осуществлялась в ходе анализа записей пациентки 1964 г. рождения, страдающей патологией cortical tremor. Данные были получены за период работы авторов в Научном центре «Юлих» в ФРГ (Forschungzentrum Juelich GmbH) во время двух сеансов съёма. Данные снимались стандартным образом (с расположением 21-го электрода по схеме «10-20%»), физиологическое состояние пациента контролировалась специалистами-медиками. Были соблюдены все основные требования к съёму ЭЭГ как с точки зрения качества полученных данных, так и в отношении здоровья пациента. Письменное согласие пациента было получено у него на весь период нахождения в Научном центре.

Для таких пациентов априорно известно о существовании патологического воздействия некоторых областей коры головного мозга на конечность, вызывающего тремор этой конечности [14]. В качестве примера ниже приведены результаты анализа записей поверхностной ЭЭГ с охваченной патологией области коры и сигнала акселерометра с треморной конечности. Длительность записей составляет около 40 секунд, частота оцифровки 200 Гц. Перед обработкой сигналы ЭЭГ и акселерометра фильтровались в диапазоне 3.5-7.5 Гц, соответствующем максимуму в спектре колебательной активности руки.

Временные ряды и периодограммы акселерометрического сигнала и сигнала ЭЭГ представлены на рис. 1.

2.2. Оценка связанности по экспериментальным данным. Оценка связанности методом причинности по Грейнджеру проводилась по двум наборам экспериментальных данных как линейным методом при различных размерностях модели Ds и Da (перебирались несколько значений вплоть до весьма больших Ds = 50, Da = 50) при лаге l = 1, так и нелинейным методом с параметрами Ds = 3, Da = 1, l = 1 и аппоксимирующей функцией g в виде кубического полинома (использовать большие значения размерности и/или степени полинома затруднительно из-за ограниченности объёма данных). Значения размерностей модели Da и Ds увеличивались до тех пор, пока не начала возрастать ошибка предсказания совместной модели ej. Результаты применения представлены на рис. 2.

Из рис. 2, в, г видно, что как линейный, так и нелинейный подходы дают возможность выявить значимое воздействие в направлении от руки к коре головного мозга. Факт наличия такого воздействия известен из исследований в области физиологии [15]. Успех линейной модели в данном случае может объясняется свойствами ЭЭГ. Как было показано в [16], сигнал ЭЭГ по свойствам близок к процессу авторегрессии-скользящего среднего высокого порядка.

Связь в обратном направлении - от головного мозга к руке, выявляется только с использованием нелинейного метода грейнджеровской причинности (рис. 2, б). Оценки, сделанные линейным методом (рис. 2, a), остаются незначимыми при переборе параметров модели: лагов l и размерностей Ds и Da в широких пределах. Между тем, наличие связи в этом направлении известно априорно из физиологических соображений [17,18]. Для объяснения данного результата можно отметить, что акселерометрический сигнал похож на сигнал автоколебательной, то есть нелинейной системы и, видимо, не может быть описан линейной моделью в достаточной степени.

Рис. 1. Временные ряды экспериментальных сигналов: а - сигнал ЭЭГ, фильтрованный в диапазоне 3.5-7.5 Гц; б - нефильтрованный акселерометрический сигнал; в, г - их периодограммы соответственно

Рис. 2. Зависимость улучшения прогноза PI, нормированного на полный 95%-й квантиль PIabs от времени задержки А (сплошная чёрная линия): a, в - для линейной модели, б, г - для нелинейной; а, б - в направлении от коры головного мозга (ЭЭГ) к руке (сигнал акселерометра) и в, г - от руки к коре. Штриховая серая линия показывает нормированный полный 95%-й квантиль. Значения, превосходящие 1, считаются значимыми с вероятностью 95%

3. Численное моделирование

3.1. Модельные системы. В ходе численного моделирования анализировались сигналы стохастических осцилляторов, связанных в общем случае с задержкой. Параметры систем выбирались с целью качественного воспроизведения свойств экспериментальных временных рядов (см. раздел 2.). Для описания колебательных свойств сигнала треморной конечности использовался осциллятор ван дер Поля

, 2\ ¿У , 2

~ " 7 = "V(Ь) +

+ уху (х(г - АХу) - у (г))3, (4)

- (т - у2) ¿г + у = "V (г) + кху (х(г - Аху) - у (г)) +

в периодическом режиме с параметрами г = 0.05 и шу = 1 с шумом (Ь) (нормальный некоррелированный шум со среднеквадратичным отклонением о,и = 0.5 и нулевым средним) (рис. 3, а, в).

Для качественного моделирования сигнала поверхностной энцефалограммы в соответствии с [16] был взят возбуждаемый шумом " (нормальный некоррелированный шум со среднеквадратичным отклонением 01 = 0.5 и нулевым средним) линейный диссипативный осциллятор с параметрами у = 0.15 и юх = 1 (рис. 3, б, г).

—х + 2у — + «Хх = ш + кух{ у (г - Аух) - х(г)) +

¿1 х ¿х 2

& +2уТг + ^ з

+ уух (у(г - Аух) - х(г)) , (5)

где Аух и Аху - временные задержки распространения сигнала между системами (4) и (5).

Уравнениями линейного осциллятора, как показано в работе [16], может быть в первом приближении описана динамика поверхностной энцефалограммы. Системы (4) и (5) решались численно методом Эйлера с шагом Н = 0.01. Для сопоставления результатов с экспериментом полагали, что Н = АЬ = 0.005 с, длина временного ряда бралась равною 8000 значений. Методы применялись при различных значениях параметров связи.

3.2. Результаты тестирования на моделях. На рис. 4 представлены результаты применения методов линейной и нелинейной грейнджеровской причинности для выявления связей между (4) и (5) для Аху = 50 мс и Аух = 0 и параметров связи Кух = 0.15, Кху = 0.05, уух = 0.05 и уух = 0.15. При этом, несмотря на то, что линейная компонента воздействия присутствует в связи в обоих направлениях, в направлении от (5) к (4) линейная модель даже достаточно высокой размерности (Д, = 50 и Оа = 50) не может выявить значимую связь; в то время как нелинейная модель низкой размерности (П,3 = 3, Оа = 1) с полиномом 3-го порядка не только обнаруживает её, но и позволяет достаточно точно определить значение задержки Аху. Аналогичные результаты были получены и при других значениях параметров

при различных уровнях шума (значениях 02 и о2и).

Рис. 3. а, б - временные ряды модельных систем (4) и (5), в, г - их спектры мощности

Рис. 4. Зависимость улучшения прогноза Р1, нормированного на полный 95%-ный квантиль Р1аЬе, от временной задержки А (чёрная сплошная линия): а, в - для линейной модели; б, г - для нелинейной; а, б - в направлении от (5) к (4); в г - от (4) к (5). Штриховая серая линия обозначает 95%-й квантиль (равен 1 вследствие нормализации). Значения, превышающие 1, признаются значимыми с вероятностью 95%

Заключение

Анализ записей ЭЭГ и сигнала акселерометра пациента показал, что воздействие конечности на кору при патологии cortical tremor выявляется как линейным, так и нелинейным методами; вместе с тем, ожидаемое из априорной информации воздействие коры головного мозга на конечность диагностируется только при использовании нелинейного подхода. Полученные оценки связей между активностью мозга и движениями руки могут указывать на принципиально нелинейный характер взаимодействия моторной коры головного мозга с конечностями. Достоверность этого вывода подтверждается не только непосредственной проверкой значимости с помощью суррогатных данных, но и хорошем соответствием с результатами оценки связанности модельной системы с нелинейной связью.

Эти выводы не контрастируют с результатами более ранних исследований (например, [19]), где связь между головным мозгом и конечностью при данной патологи обнаруживалась линейными методами, например, с помощью анализа когерентности. Однако использованные подходы не позволили сделать выводы о её направленности. Показательно, что этот вывод соответствует существующим представлениям физиологов о механизме рассматриваемой патологии.

Таким образом, результаты данной работы подтверждают целесообразность применения нелинейных моделей в задаче оценки связанности между системами по их экспериментальным сигналам несмотря на дополнительные проблемы, возникающие при использовании нелинейного подхода, такие как зависимость результатов от конкретного вида нелинейности и чрезмерное увеличение числа коэффициентов при использовании аппроксимирующих функций общего вида (алгебраических и тригонометрических полиномов), что, в свою очередь, повышает влияние шумов и требования к длине временной реализации.

Работа поддержана Российским фондом фундаментальных исследований -грант 08-02-00081 и Аналитической ведомственной целевой программы «Развитие научного потенциала высшей школы (2009-2010 годы)» - проект 2.1.1/1738.

Библиографический список

1. Флейшман А.Н.Вариабельность ритма сердца и медленные колебания гемодинамики (нелинейные феномены в клинической практике). Новосибирск: Изд-во Сибирского отделения РАН, 2009. 193 с.

2. Baccala L.A., Sameshima K., Ballester G., Do Valle A.C. and Timo-Laria C. Studing the interactions betweenbrain structures via directed coherence and Granger causality // Applied Sig. Processing. 1998. Vol. 5. P. 40.

3. Nalatore H., Ding M., Rangarajan G. Mitigating the effects of measurement noise on Granger causality // Phys. Rev. E. 2007. Vol. 75. 031123

4. Feldman U. and Bhattacharya /.Predictability improvement as an asymmetrical measure of interdependence in bivariate time series // Int. J. of Bifurcation and Chaos. 2004. Vol. 14, № 2. P. 505.

5. Granger C.W.J. Investigating causal relations by econometric models and cross-spectral methods // Econometrica. 1969. Vol. 37, № 3. P. 424.

6. 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.

7. Marinazzo D., Pellicoro M., and Stramaglia S. Nonlinear parametric model for Granger causality of time series // Phys. Rev. E. 2006. Vol. 73. 066216.

8. Ishiguro K., Otsu N., Lungarella M. and Kuniyoshi Y Comparison of nonlinear Granger causality extensions for low-dimensional systems // Phys. Rev. E. 2008. Vol. 77. 036217.

9. Rosenblum M.G. and Pikovsky A.S. Detecting direction of coupling in interacting oscillators // Phys. Rev. E. 2001. Vol. 64. 045202.

10. Smirnov D. and Bezruchko B. Estimation of interaction strength and direction from short and noisy time series // Phys. Rev. E. 2003. Vol. 68. 046209.

11. SmirnovD., Barnikol U.B., Barnikol T.T., Bezruchko B.P., Hauptmann C., Buehrle C., Maarouf M., Sturm V., Freund H.-J., and Tass P.A. The generation of parkinsonian tremor as revealed by directional coupling analysis // Europhysics Letters. 2008. Vol. 83. 20003.

12. Dolan K. and Neiman A. Surrogate analysis of coherent multichannel data // Phys. Rev. E. 2002. Vol. 65. 026108.

13. Schreiber T. and Schmitz A. Surrogate time series // Physica D. 2000. Vol. 142. 346.

14. Hallett M. Overview of human tremor physiology // Movement Disorders. Vol. 13, Issue S3. P. 43

15. Wang S., Aziz T., Stein J., Liu X. Time-frequency analysis of transient neuromuscular events: dynamic changes in activity of the subthalamic nucleus and forearm muscles related to the inter-mittent resting tremor // J. Neurosci. Methods. 2005. Vol. 145. P. 151.

16. Timmer J., Haussler S., Lauk M., Lucking C.-H. Pathological tremor: deterministic chaos or nonlinear stochastic oscillators // Chaos. 2000. Vol. 10, № 1. P. 278.

17. Eichler M. Graphical modelling of dynamic relationships in multivariate time series // Handbook of Time Series Analysis / Ed. M. Winterhalder, B. Schelter and J. Timmer. Berlin: Wiley-VCH, 2006. P. 335.

18. Toro C., Pascual-Leone A., Deuschl G., Tate E., Pranzatelli M.R., and Hallett M. A common manifestation of cortical myoclonus // Neurology. 1993. Vol. 43(11). P. 2346.

19. Hellwig B., Haufiler S., Schelter B., Lauk M., Guschlbauer B., Timmer J., Lucking C.H. Tremor-correlated cortical activity in essential tremor // The Lancet. 2001. Vol 357. P. 519.

СФ Института радиотехники Поступила в редакцию 16.07.2010

и электроники им. Котельникова РАН После доработки 21.09.2010

ROLE OF MODEL NONLINEARITY FOR GRANGER CAUSALITY BASED COUPLING ESTIMATION FOR PATHOLOGICAL TREMOR

I.V. Sysoev, A.S. Karavaev, P.I. Nakonechny

Estimating coupling between systems of different nature is an urgent field of nonlinear dynamics method application. This work aims to compare classical linear Granger

approach and its nonlinear analogues based on analysis of ethalon dynamical systems and neurophysiological data. The results achieved show nonlinear approach to be more sensitive, and so it is able to detect significant coupling, when linear one fails.

Keywords: Coupling analysis, nonlinear Granger causality, processing EEGs, dynamic modeling.

Сысоев Илья Вячеславович - родился в Саратове (1983). Окончил Лицей прикладных наук (1999) и Факультет нелинейных процессов (2004). Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2007). Работал на кафедре электроники, колебаний и волн (2005-2007), в настоящее время доцент базовой кафедры динамического моделирования и биомедицинской инженерии. Научные интересы - исследование сигналов биологической природы методами нелинейной динамики, исследование эффективности и модернизация подходов к анализу сигналов. Автор более 40 публикаций.

410012, Саратов, ул. Астраханская, 83

Саратовский государственный университет им. Н.Г. Чернышевского E-mail: ivssci@gmail.com

Караваев Анатолий Сергеевич - родился в 1981 году. Окончил факультет нелинейных процессов СГУ (2004). Кандидат физико-математических наук (2007). Доцент базовой кафедры динамического моделирования и биомедицинской инженерии ФНБМТ СГУ, с.н.с. СФ ИРЭ РАН. Научные интересы -экспериментальное исследование нелинейных явлений в системах различной физической природы. Автор около 20 статей в отечественных и зарубежных научных журналах.

410012, Саратов, ул. Астраханская, 83

Саратовский государственный университет им. Н.Г. Чернышевского E-mail: karavaevas@gmail.com

Наконечный Павел Игоревич - родился в Саратове (1985). Окончил Физико-технический лицей (2002) и факультет нелинейных процессов СГУ (2007). В настоящее время аспирант кафедры динамического моделирования и биомедицинской инженерии. Научные интересы - анализ временных рядов, информационные технологии.

410012, Саратов, ул. Астраханская, 83

Саратовский государственный университет им. Н.Г. Чернышевского E-mail: zorg1331@gmail.com