CC BY
42
Изв. вузов «ПНД», т. 15, № 4, 2007 УДК 530.182(075.8)
ПРИМЕНЕНИЕ НЕПРЕРЫВНОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ ДЛЯ АНАЛИЗА ПЕРЕМЕЖАЮЩЕГОСЯ ПОВЕДЕНИЯ
А.А. Короновский, И.М. Минюхин, А.А. Тыщенко, А.Е. Храмов, И.С. Мидзяновская, Е.Ю. Ситникова, G. van Luijtelaar, C.M. van Rijn
В данной работе предлагается эффективный метод анализа сигналов при помощи непрерывного вейвлет-преобразования. Рассматривается применение данного метода для определения длительности ламинарных и турбулентных фаз движения для перемежающегося поведения различных типов, включая анализ временных рядов, порожденных живыми системами. Показано, что предложенный метод обладает высокой устойчивостью к шумам и флуктуациям, искажающим исходную временную реализацию.
Введение. Обсуждение проблемы
Многие процессы в природе и физике носят перемежающийся характер. В частности, в гидродинамике в ряде случаев при больших числах Рейнольдса можно наблюдать перемежающуюся структуру течения [1]. Перемежаемость является одним из классических сценариев перехода от периодических колебаний к хаотическим [1,2]. При этом сигнал представляет собой последовательность чередующихся регулярных (ламинарных) фаз и нерегулярных всплесков (турбулентные фазы). Число нерегулярных всплесков нарастает при увеличении управляющего параметра до тех пор, пока движение в системе полностью не хаотизируется. Перемежающееся поведение наблюдается вблизи границы возникновения режимов хаотической синхронизации связанных осцилляторов [3-5]; чередование судорожной активности и «нормального» функционирования мозга у животных, генетически предрасположенных к абсанс эпилепсии, также представляет собой перемежаемость [6]. Существует определенная классификация перемежающегося поведения, в частности, выделяют перемежаемость типа I-III [7, 8], on-off перемежаемость [9], перемежаемость типа «игольное ушко» [10].
При исследовании перемежаемости важной проблемой является задача о выделении в сигналах временных интервалов, соответствующих различным типам динамики систем, демонстрирующих перемежающееся поведение (задача о диагностике ламинарных и турбулентных фаз). Как правило, участки временной реализации, на которых наблюдается поведение, близкое к регулярному (так называемые ламинарные фазы), выделяются легче, нежели участки с нерегулярным (турбулентным)
поведением. Классическими методами выделения регулярных фаз движения являются методы, основанные на анализе «текущего» периода колебаний или амплитуды системы. Первый метод может работать лишь тогда, когда ламинарная фаза представляет сигнал, очень близкий к строго регулярному, что имеет место далеко не всегда. Второй метод может быть применен только в том случае, когда турбулентная фаза характеризуется существенно большей амплитудой колебаний по сравнению с ламинарной фазой. В этом случае обычно поступают следующим образом: для наблюдаемой величины задается определенное пороговое значение А и пока состояние системы характеризуется значением, по абсолютной величине превышающим этот порог, считается, что в системе реализуется турбулентная фаза, а когда значение, описывающее состояние системы, оказывается ниже порогового - ламинарная. Подобный подход применялся, например, при анализе перемежающегося поведения на границе синхронизации с запаздыванием (см., например, [11]). В этом случае анализируемый сигнал и(Ь) представлял собой разность состояний двух взаимодействующих систем, сдвинутых друг относительно друга на время запаздывания т: и(Ь) = Х\(Ь) — Х2(Ь — т). В те моменты времени, когда поведение систем было синхронизовано, величина и(Ь) была близка к нулю (ламинарная фаза), а когда синхронизация разрушалась (турбулентная фаза), величина и(Ь) изменялась в достаточно широких пределах. Однако в данном случае недостаток метода заключается в том, что в течение турбулентной фазы движения величина и(Ь) может пересекать нулевое значение, то есть оказываться как выше, так и ниже порога А. В свою очередь, это приводит к погрешности выделения турбулентных фаз, так как каждая турбулентная фаза оказывается разделенной на несколько частей и, плюс ко всему, появляются ошибочно диагностированные короткие ламинарные фазы. Таким образом, получается погрешность при анализе статистических свойств перемежающегося поведения (определение средней длительности, законов распределения ламинарных и турбулентных фаз).
В работе [12] для исследования перемежающегося поведения предлагалось использовать непрерывное вейвлет-преобразование, которое является эффективным средством анализа временных реализаций, характеристики которых изменяются с течением времени [13-15]. Метод, предложенный в [12], хорошо работает только если поведение системы в ламинарной фазе является регулярным и имеет одну и ту же структуру в любой момент времени ламинарной фазы системы. Суть метода заключается в следующем: определяется число максимумов вейвлет-поверхности в каждый момент времени Ь, и интервалы времени, в течение которых число максимумов на вейвлет-поверхности не изменяется, считаются ламинарными фазами, а другие, соответственно, - турбулентными. Однако данный метод оказывается неприменимым для систем, у которых амплитуда сигнала изменяется в широком диапазоне, а ламинарная фаза непериодична (в этом случае число максимумов в течение ламинарной фазы для любого момента времени различно). В частности, такой режим наблюдается при анализе перехода к синхронному поведению хаотических осцилляторов [11]. Для данного класса связанных хаотических систем также был предложен метод анализа [16], однако и он имел жесткие ограничения - для его применения было необходимо, чтобы во время турбулентной фазы (асинхронная динамика связанных осцилляторов) основная энергия вейвлет-спектра Е(в, Ь) колебаний приходилась на два основных характерных временных масштаба в 1,2, что также оказалось неприменимым для широкого круга систем.
В настоящей работе предложен новый универсальный метод выделения ламинарных и турбулентных фаз, основанный на непрерывном вейвлет-преобразовании и позволяющий проводить эффективный анализ временных реализаций систем, демонстрирующих явление перемежаемости. Следует отметить, что данный метод может быть с успехом применен для анализа различных типов перемежающегося поведения в разных системах (как с непрерывным, так и с дискретным временем), включая живые системы.
Структура настоящей работы следующая: в разделе 1 описывается предлагаемый метод. Применение этого метода для анализа перемежающегося поведения различных систем рассматривается в разделе 2: перемежаемость I типа в динамической системе с непрерывным временем (система Лоренца) рассматривается в подразделе 2.1, а тот же самый тип перемежаемости в системе с дискретным временем (логистическое отображение) - в подразделе 2.2. Явление оп-оГГ перемежаемости, наблюдаемое вблизи границы возникновения режима синхронизации с запаздыванием, проанализировано с помощью предложенного метода в подразделе 2.3. В разделе 3 рассматривается устойчивость предложенного метода по отношению к шумам и флуктуациям, существующим в анализируемых сигналах. Наконец, использование предлагаемого метода для анализа перемежающегося поведения в спонтанной неконвульсивной судорожной активности у крыс изложено в разделе 4. Обобщающие замечания и выводы приведены в Заключении.
1. Метод выделения ламинарных и турбулентных фаз
Предлагаемый в работе метод выделения различных фаз в поведении системы, демонстрирующей перемежающееся поведение, основан на непрерывном вейвлет-преобразовании анализируемого временного сигнала х(1) [13]
w(1, в) = -1=1 хр)г(^^ (1)
где 'ф(п) - базисный вейвлет (звездочка обозначает комплексное сопряжение), в - временной масштаб. В качестве базисного использовался вейвлет Морле
Ц(ц) = к—1/4е^°пе—п2/2. (2)
Параметр вейвлет-функции юо был выбран равным 2п, что, с одной стороны, обеспечивает хорошее соотношение между локализациями вейвлет-функции во времени и фурье-пространстве [13], а с другой стороны, позволяет легко сопоставлять временные масштабы в вейвлет-преобразования (1) с частотами / спектрального представления сигнала, поскольку для данного значения параметра юо выполняется соотношение в & 1// [13,17].
По аналогии со спектром мощности фурье-преобразования можно ввести в рассмотрение мгновенное
Е (в,1) = \W (в,1)\2 (3)
и интегральное
т
Е (в)) = \ш (в,ь)\Чг (4)
о
распределение энергии по временным масштабам [13].
Поскольку поведение системы во время турбулентных и ламинарных фаз различается1 , то и структура вейвлет-поверхности Ш(Ь, в) в области ламинарных и турбулентных фаз движения также будет существенно различна [12, 13]. Иными словами, энергия вейвлет-спектра Е(в,Ь) будет распределена по характерным временным масштабам в, которые для разных фаз временной реализации х(Ь) будут разными, причем доля энергии, приходящейся на эти характерные временные масштабы, также будет различаться. Таким образом, можно перейти от анализа структуры вейвлет-поверхности Ш(в, Ь) к анализу распределения энергии вейвлет-спектра по характерным временным масштабам. Для выделения ламинарных и турбулентных фаз в каждый момент времени Ь определяется суммарное значение энергии вейвлет-спектра <м(Ь), приходящейся на выбранный диапазон характерных временных масштабов в е Б = (вх; в2)
-ш(г) = ! Е(Ь, в) (1в. (5)
Диапазон характерных временных масштабов в, по которым будет находиться величина <м(Ь), определяется рассматриваемой системой, и в каждом конкретном случае должен выбираться на основании мгновенных распределений энергии вейвлет-спектра для различных фаз перемежающегося поведения. В ряде случаев может сложиться ситуация, когда необходимо рассматривать несколько диапазонов временных масштабов Бг, которые однозначно позволяли бы охарактеризовать поведение системы. В этом случае интегрирование (5) должно проводиться по объединению Б = У Si соответствующих диапазонов временных масштабов.
г
2. Выделение ламинарных и турбулентных фаз в поведении нелинейных систем
Рассмотрим теперь применение предложенного метода для анализа различных типов перемежающегося поведения в эталонных нелинейных системах, таких как система Лоренца, логистическое отображение и связанные системы Ресслера.
2.1. Перемежаемость I типа в системе Лоренца. Начнем рассмотрение использования вышеописанного метода применительно к системе Лоренца [12]
Х = а (у — х),
у = гх — хг — у, (6)
г = ху — Ьг.
1Следует обратить внимание, что для различных типов перемежающегося поведения под «ламинарными» и «турбулентными» фазами будут пониматься разные типы поведения.
Известно, что при значениях управляющих параметров о = 10,6 = 8/3 и при г < г* = 166.07 временная реализация х(1) представляет собой устойчивое периодическое движение. При превышении порога г* регулярные колебания х(1) (ламинарная фаза) прерываются хаотическими всплесками, которые с ростом г становятся все более и более длительными, пока движение полностью не хаотизируется. Перемежаемость в данной системе классифицируется как перемежаемость I типа [18].
На рис. 1 показаны временная реализация х(1) системы Лоренца (6) при значении бифуркационного параметра г = 166.14, соответствующая ей проекция вейвлет-поверхности \W(1, в)\ на плоскость (в,1) и энергия вейвлет-спектра ,ш(1), приходящаяся на диапазон характерных временных масштабов в Е (в1; в2) (в1 = 0.15, в2 = 2.4). Наиболее темные участки проекции \ W(1, в)\ на плоскость (в, 1) (рис. 1, б) соответствуют максимумам вейвлет-поверхности, кривая на рисунке ограничивает область влияния краевых эффектов. На проекции вейвлет-поверхности четко выделяются структуры, соответствующие ламинарной и турбулентной фазам перемежающегося поведения во временной реализации. Структура вейвлет-поверхности, соответствующая регулярной фазе движения, имеет профиль W(1, в)\г=соп^ с двумя глобальными максимумами, соответствующими характерным временным масштабам ва & 0.4 и вь & 1.1 (см. рис. 2, а, на котором показано мгновенное распределение энергии вейвлет-спектра по временным масштабам Е(в)). Этим максимумам вейвлет-поверхности соответствуют две темные области на проекции поверхности
Рис. 1. а - Временная реализация х(£) системы Лоренца (6) при г = 166.14. б - Проекция вейвлет-поверхности (£,в)\ на плоскость (£,в) для временной реализации х(£). Интенсивность градаций серого цвета пропорциональна абсолютной величине вейвлет-поверхности (в,£)|. Сплошная линия ограничивает область влияния краевых эффектов, достоверные данные располагаются ниже этой линии [13, 15]. в - Энергия вейвлет-спектра ш(£), приходящаяся на диапазон характерных временных масштабов в е (0.15, 2.4)
\W(t, s)\, расположенные параллельно оси времени (см. рис. 1, б). Профиль поверхности \W(t, s)\ не изменяется с течением времени t в пределах ламинарной фазы. Со входом системы в турбулентную фазу вид поверхности \ W(t, s)\ сильно изменяется. Можно сказать, что в режиме хаотической динамики наблюдается «всплеск» разномасштабных колебательных явлений. Это хорошо видно на рис. 2, б: распределение энергии вейвлет-спектра по временным масштабам E(s) существенным образом отличается от распределения, соответствующего ламинарной фазе (ср. с рис. 2, a). Из рис. 1, б видно, что области вейвлет-поверхности, соответствующие турбулентным фазам, четко локализованы во времени [12]. Также видно, что в пределах ламинарной фазы функция w(t) имеет структуру, близкую к периодической, и совершает колебания малой амплитуды. В течение турбулентной фазы w(t) испытывает резкий провал (см. рис. 1, в), что позволяет точно диагностировать различные фазы движения. Моменты времени, в которые величина w(t) оказывается больше порогового значения А = 250, считаются ламинарными фазами, а интервалы времени, для которых выполняется условие w(t) < А - турбулентными.
На рис. 3 представлена зависимость средней длительности ламинарной tb и турбулентной tx фаз от параметра надкритичности (r — r*). Благодаря логарифми-
Рис. 2. Профили вейвлет-поверхности \W(s,t)\, соответствующие ламинарной (a) и турбулентной (б) фазам временной реализации x(t) системы Лоренца. На профиле вейвлет-поверхности ламинарной фазы отчетливо видны два максимума, соответствующие характерным временным масштабам sa ~ 0.4 И Sf, ~ 1.1
îlJT 8.0
6.0 4.0 H 2.0 0.0
h
□ □ □ □ □
0.1 0.2 0.4 0.S
Рис. 3. а - Зависимость средней длительности ламинарных фаз от параметра надкритичности вблизи порога перемежаемости для системы Лоренца (6); пунктирной линией показана степенная зависимость с показателем степени а = -1/2. б - Средние длительности ламинарной (•) и ¿т турбулентной (□) фаз колебаний х(Ь) в системе Лоренца вблизи порога перемежаемости в зависимости от параметра надкритичности
ческому масштабу (рис. 3, а) видно, что средняя длительность ламинарной фазы зависит от параметра надкритичности как tb ~ (r — r*)-1/2. На рис. 3, б показаны зависимости средней длительности ламинарной (•) и турбулентной (□) фаз. Видно, что средняя длительность хаотической фазы практически не меняется с увеличением значения параметра надкритичности. Это объясняется тем, что турбулентная фаза движения есть не что иное, как стадия реламинаризации, то есть возвращение к регулярному (ламинарному) движению. Понятно, что длительность реламинариза-ции определяется системой и практически не зависит от параметра надкритичности. Следует отметить, что полученные результаты по зависимости средней длительности ламинарной фазы от параметра надкритичности хорошо согласуются с теорией перемежаемости I типа [18]. Следует также отметить, что соответствующая зависимость средней длительности турбулентных фаз от параметра надкритичности, как правило, не анализируется, что связано как со сложностью получения аналитических результатов, так и с трудностями выделения турбулентных фаз из анализируемой временной реализации.
2.2. Перемежаемость I типа в логистическом отображении. Рассмотрим теперь применение предложенного метода для выделения ламинарных и турбулентных фаз перемежающегося поведения I типа в системе с дискретным временем. Для этого рассмотрим логистическое отображение вида
Xn+i = exn(1 — xn). (7)
При значении управляющего параметра e* = 1 + V8 в системе происходит касательная бифуркация цикла периода 3, а ниже по параметру ламинарные фазы (колебания, близкие к периодическим с периодом три) прерываются турбулентными (хаотическими) всплесками.
На рис. 4 показана временная реализация xn логистического отображения для значения параметра надкритичности (e* — e) = 0.001, соответствующая ей проекция вейвлет-поверхности \W(n,s) \ и функция
Wn = У \W(s,n)\2 ds, (8)
S
построенная для диапазона временных масштабов s Е (s1; s2) (s1 = 2, s2 = 4). Из рисунка видно, что в данном случае структура вейвлет-поверхности в течение ламинарных и турбулентных фаз движения, так же как и для системы Лоренца (6), существенно различается, и функция Wn позволяет достаточно точно диагностировать ламинарные и турбулентные фазы, точно так же, как и для системы с потоковым временем (см. подраздел 2.1.).
Понятно, что все закономерности, характерные для перемежаемости I типа, будут проявляться и для системы с дискретным временем. На рис. 5 представлена зависимость средней длительности ламинарной фазы от параметра надкритичности (e—e*), полученная для логистического отображения. Зависимость средней длительности ламинарной фазы от параметра надкритичности также подчиняется степенному закону с показателем степени а = —1/2, что полностью соответствует теории перемежаемости I типа [18]. Средняя длительность турбулентной фазы практически не меняется с ростом параметра надкритичности (рис. 5, б).
Рис. 4. а - Временная реализация xn логистического отображения (7) при значении параметра надкритичности (s* — s) = 0.001 и б - проекция вейвлет-поверхности \W(n,s)\ на плоскость (n,s) для временной реализации xn. Интенсивность градаций серого цвета пропорциональна абсолютной величине вейвлет-поверхности \W(s, t)\. Сплошная линия ограничивает область влияния краевых эффектов, достоверные данные располагаются ниже этой линии. в - Мгновенная энергия вейвлет-спектра wn, приходящаяся на диапазон характерных временных масштабов s G (2.0, 4.0) в моменты дискретного времени n, пороговое значение было выбрано А = 0.4
Рис. 5. a - Зависимость средней длительности ламинарных фаз от параметра надкритичности вблизи порога перемежаемости для логистического отображения (7). Пунктирная линия соответствует степенной функции с показателем а = —1/2. б - Средние длительности tL ламинарной (•) и tT турбулентной (□) фаз колебаний xn для логистического отображения вблизи порога перемежаемости в зависимости от параметра надкритичности (s — s*)
2.3. Перемежающаяся синхронизация с запаздыванием в системе связанных осцилляторов Ресслера. До этого момента рассматривались системы (как с непрерывным, так и дискретным временем), в которых наблюдался режим перемежаемости I типа. Покажем теперь применимость данного метода для систем, в которых реализуется перемежающееся поведение принципиально другого класса. Рассмотрим систему осцилляторов Ресслера, связанных двунаправленной связью [3]
х 1,2 = -Ы1,2У1,2 - ¿1,2 + е(х2,1 - х1,2),
У/1,2 = ®1,2х1,2 + ау1,2, (9)
-¿1,2 = f + ¿1,2(х1,2 - С),
где ю1 = 0.99, ю2 = 0.95, а = 0.165, f = 0.2, с = 10. При увеличении параметра связи е в системе поочередно реализуются режимы фазовой синхронизации, синхронизации с запаздыванием, и при дальнейшем увеличении параметра связи е системы стремятся к режиму полной синхронизации [19,20].
Режим полной синхронизации характеризуется выполнением для векторов состояния взаимодействующих систем соотношения \х1 (1) — Х2 (1)\ & 0. Режим синхронизации с запаздыванием [4] характеризуется тем, что колебания в одной системе отстают от колебаний в другой на некоторый временной сдвиг т: Х2(1 + т) & Х1(1). С увеличением параметра связи между хаотическими системами временной сдвиг т уменьшается и стремится к нулю, а системы стремятся к режиму полной синхронизации. Переход к режиму синхронизации с запаздыванием сопровождается режимом перемежающейся синхронизации с запаздыванием. Для рассматриваемых систем (9) с указанным набором значений управляющих параметров перемежающееся поведение наблюдается при параметре связи 0.1 <е< 0.15 [4, 11]. В этом случае имеет смысл рассматривать в качестве исследуемого сигнала разность состояний и(1) = Х2(1 + т)х1(1), абсолютная величина которой в случае ^-синхронизации стремится к нулю (\и(1)\ ^ 0). В режиме перемежающейся синхронизации с запаздыванием данный сигнал представляет собой последовательность чередующихся участков с малой амплитудой (ламинарных фаз), соответствующих интервалам времени синхронной динамики связанных систем, и нерегулярных всплесков (турбулентных фаз), которые отвечают моментам времени, когда системы демонстрируют асинхронное поведение. Число нерегулярных всплесков уменьшается с увеличением параметра связи. Перемежаемость в данном случае классифицируется как вп-в// перемежаемость [9,21].
Для определения времени запаздывания т обычно используют функцию подобия (см. например [3,4,11]). Это не совсем удобно, поскольку для различных значений параметра связи е величины запаздывания т будут различаться, и, соответственно, процедуру нахождения т с помощью функции подобия необходимо выполнять для каждого значения параметра связи. В работах [22,23] было показано, что величина т для данной системы пропорциональна сдвигу между основными спектральными компонентами спектра Фурье и в зависимости от параметра связи е подчиняется степенному закону с показателем а = —1. Для системы Ресслера (9) с выбранными значениями управляющих параметров было получено т = 0.0418е_1 [16].
На рис. 6 показаны временная реализация и(1) = х2(1 + т) — х1(1) при значении параметра связи е = 0.1, соответствующая ей проекция вейвлет-поверхности
Рис. 6. а - Временная реализация и(Ь) = х2(Ь + т) — х\(Ь) системы двух связанных осцилляторов Рес-слера (9) при значении параметра связи е = 0.1 и б - проекция вейвлет-поверхности (Ь,в)\ на плоскость (Ь, в) для временной реализации и(Ь). По оси абсцисс отложено время Ь, по оси ординат -временные масштабы в. Интенсивность градаций серого цвета пропорциональна абсолютной величине вейвлет-поверхности (в,Ь)\. Сплошная линия ограничивает область влияния краевых эффектов, достоверные данные располагаются ниже этой линии. в - Энергия вейвлет-спектра 'ш(Ь), приходящаяся на диапазон характерных временных масштабов в е (5,15)
(Ь, в)\ и функция w(t), характеризующая значение мгновенной энергии вейвлет-спектра, приходящейся на диапазон характерных временных масштабов в Е (81; в2), в! = 5, в2 = 15. Вводя в рассмотрение некоторое пороговое значение А, можно легко выделить ламинарные и турбулентные фазы в перемежающемся поведении системы (9), которые соответствуют синхронной и асинхронной динамике связанных систем Ресслера в режиме перемежающейся синхронизации с запаздыванием. При анализе функции w(t) были получены зависимости средней длительности ламинарных и турбулентных фаз от параметра надкритичности (е* — е), где е* ~ 0.15 -порог установления режима синхронизации с запаздыванием. При увеличении параметра надкритичности в диапазоне 0.1 < е < 0.15 число турбулентных (а соответственно, и ламинарных) фаз резко уменьшается, ламинарные фазы становятся более длительными, а турбулентные - более короткими. На рис. 7, а в двойном логарифмическом масштабе показана зависимость средней длительности Ьь ламинарных фаз от величины параметра надкритичности (е* — е). Видно, что средняя длительность ламинарных фаз подчиняется степенному закону Ьь ~ (е* — е)_1, что характерно для режима оп-оГГ перемежаемости. На рис. 7, б приведена аналогичная зависимость средней длительности Ьт турбулентных фаз от параметра надкритичности (е* — е) (по оси ординат - масштаб логарифмический). Из рисунка следует, что средняя дли-
5000 1000 200 ■
0.005 0.01 0.02 £*-£ _ 0.005 0.025 £*-£
а о
Рис. 7. а - Зависимость средней длительности ламинарных фаз от параметра надкритичности вблизи порога возникновения режима синхронизации с запаздыванием для системы (9). Пунктирная линия -степенная функция с показателем а = —1. Зависимость приведена в двойном логарифмическом масштабе. б - Зависимость средней длительности турбулентной фазы от параметра надкритичности. Пунктирная линия - экспоненциальная функция с показателем, зависящим от параметра надкритично-сти (е* — е). Масштаб по оси абсцисс - логарифмический, по оси ординат - линейный
тельность 1т турбулентных фаз зависит от параметра надкритичности (е* — е) по экспоненциальному закону. Следует отметить, что данные результаты хорошо согласуются с полученными ранее [16].
3. Влияние флуктуаций и шумов
Как видно из раздела 2, предложенный метод позволяет очень хорошо выделять ламинарные и турбулентные фазы при различных типах перемежающегося поведения эталонных систем. В то же самое время, почти все временные ряды, полученные в ходе экспериментальных измерений, искажаются за счет влияния различных факторов. Поэтому чрезвычайно важной проблемой, тесно связанной с выделением фаз различной динамики в перемежающемся поведении систем, является учет влияния шумов. Очевидно, что воздействие шумов может приводить к ошибкам в выделении ламинарных и турбулентных фаз, причем число таких ошибок должно расти с увеличением интенсивности шума. В связи с этим возникает вопрос об устойчивости предложенного метода выделения различных фаз перемежающегося поведения к влиянию шумов и флуктуаций. Рассмотрим данный вопрос более подробно.
В качестве исследуемой модели будем рассматривать систему Лоренца (6) со значениями управляющих параметров, указанных в подразделе 2.1. К временной реализации системы Лоренца х(1) будем добавлять случайную величину ^(1) с нулевым средним, равномерно распределенную на интервале (—0.5; 0.5). Численное моделирование случайной величины ^(1) осуществлялось так, как это описано в [24]. Интенсивность шума будем характеризовать величиной В
хп (1)= х(1) + В£(1). (10)
Случай В = 0 соответствует поведению системы без шума, выделение ламинарных и турбулентных фаз для которого было подробно рассмотрено в подразделе 2.1. Очевидно, что увеличение В будет приводить к искажениям исходной временной реализации, и, соответственно, ламинарные и турбулентные фазы, определенные
по такому зашумленному временному ряду хв (Ь), будут отличаться от фаз, выделенных по временной реализации х(Ь) исходной системы (6) без шума. Соответственно, исходная временная реализация х(Ь) может рассматриваться как «эталонная», и, сопоставляя ламинарные и турбулентные фазы, выделенные по временным реализациям х(Ь) и хв(Ь), можно сделать вывод о том, к каким ошибкам в выделении турбулентных и ламинарных фаз приводит добавление шума.
Отметим, что искажения могут быть двух видов. Во-первых, в зашумленной временной реализации хв (Ь) может быть обнаружена турбулентная фаза, которой не существует во временной реализации без шума, или наоборот, не найдена турбулентная фаза, присутствующая в исходной временной реализации х(Ь). Во-вторых, турбулентная фаза по зашумленному временному ряду хв (Ь) может быть выделена, но ее длительность, моменты времени возникновения и окончания будут отличаться от истинных, определенных по временной реализации системы без шума х(Ь).
На рис. 8 показаны временные реализации х(Ь) и хв(Ь) системы Лоренца (6) без добавления шума (рис. 8, а) и при наличии шума с интенсивностью О = 48 (рис. 8, б), соответственно. Из рисунка видно, что добавление шума серьезным об-
Рис. 8. Временные реализации системы Лоренца (6) для значения управляющего параметра г = 166.14: а - х(Ь) без шума (Б = 0); б - хв (Ь) с добавлением шума (10), интенсивность которого Б = 48 соответствует примерно половине амплитуды исходного временного ряда. Турбулентные фазы, выделенные по данным временным реализациям методом, описанным в разделе 1, показаны квадратными скобками: пунктирные скобки соответствуют турбулентным фазам, выделенным по временной реализации без шума х(Ь), сплошные - по временной реализации хв(Ь) с шумом интенсивности Б = 48
разом искажает временной ряд, так что визуально выделить ламинарные и турбулентные фазы не представляется возможным. В то же самое время, используя метод, предложенный в разделе 1, удается легко выделить ламинарные и турбулентные фазы как в исходной, так и в зашумленной временной реализации. На рис. 8 выделенные турбулентные фазы показаны скобками. Из рисунка видно, что, несмотря на значительную интенсивность добавленного шума В^(1), предложенный метод позволяет весьма эффективно выделить ламинарные и турбулентные фазы, причем влияние шума приводит, прежде всего, к незначительному изменению диагностирования моментов начала и конца ламинарной фазы. Следует также отметить, что при выделении ламинарных и турбулентных фаз по временной реализации с шумом хп (1) при больших интенсивностях шума В возникают ошибочно определенные турбулентные фазы с очень короткой длительностью. Эти ошибочно найденные фазы легко могут быть исключены из рассмотрения с помощью анализа длительности турбулентных фаз2.
Для того чтобы объективно охарактеризовать ошибки, возникающие при наличии шума, введем в рассмотрение вспомогательную временную реализацию в(1), построенную по анализируемому временному ряду х(1):
{1, для турбулентной фазы,
0 " ф (11)
0, для ламинарной фазы.
Очевидно, что произведение в(1)вп (1) будет равно единице только в те моменты времени 1, когда турбулентная фаза выделена как во временной реализации хп (1) с шумом, так и в исходной временной реализации х(1) без шума. Вводя в рассмотрение оператор £[•]
+те
Ь[в(1)] = У в(1) й1, (12)
—те
можно получить суммарную длину интервалов времени, на которых анализируемая функция в(1) равна единице. Иными словами, Ь[в(1)] дает суммарную длину всех турбулентных фаз во временной реализации х(1). Очевидно, что Ь[в(1)вп(1)] даст суммарную длину всех интервалов времени, для которых турбулентная фаза была выделена как по исходной х(1), так и по зашумленной хп (1) временным реализациям. Все возникающие ошибки в выделении турбулентных фаз будут приводить к уменьшению величины Ь[в(1)вп(1)]. Тогда величина
0 Ь[з(1)зп (1)] (13)
= ьт] (13)
будет характеризовать нормированную ошибку в выделении турбулентных фаз. Понятно, что при нулевом уровне шума В = 0 величина От будет равна единице и по мере возрастания интенсивности шума В эта величина будет уменьшаться.
2В подразделе 2.1 было показано, что средняя длительность турбулентной фазы практически не зависит от величины параметра надкритичности, что определяется механизмами реламинаризации. Очевидно, что существует некоторое характерное время реламинаризации, и турбулентные фазы с очень короткой длительностью будут явным образом отличаться от всех остальных турбулентных фаз.
Аналогично, для того чтобы охарактеризовать ошибку в выделении ламинарных фаз, можно рассмотреть величину
□т
Ь[вп (1)(1 - 8(1))] Ь[1 - 8(1)] ■
(14)
Эта величина представляет собой суммарную длительность интервалов времени, на которых ламинарная фаза по временной реализации хп (1) в присутствии шума была ошибочно определена как турбулентная, нормированная на суммарную длительность ламинарных фаз, выделенных по исходной реализации х(1) без шума. При отсутствии ошибок в выделении ламинарных фаз величина Оь будет равна нулю (например, при В = 0) и по мере увеличения интенсивности шума и возрастания ошибок в определении ламинарных фаз она будет увеличиваться.
Зависимости нормированных ошибок выделения ламинарных и турбулентных фаз Оь и От приведены на рис. 9. Отчетливо видно, что с увеличением интенсивности шума В количество ошибок в определении ламинарных и турбулентных фаз увеличивается, однако заметные искажения появляются при достаточно больших интенсивностях шума, когда амплитуда случайного сигнала В'%(1) становится сопоставимой с амплитудой исходного временного ряда х(1).
Таким образом, предложенный метод обладает хорошей устойчивостью по отношению к шумам и флуктуациям и позволяет хорошо выделять ламинарные и турбулентные фазы во временной реализации системы, демонстрирующей перемежающееся поведение даже при достаточно больших интенсивностях шума. При этом основные ошибки связаны с неточным определением моментов времени начала и конца турбулентной фазы.
Рис. 9. Зависимость нормированных ошибок в определении ламинарных (пунктирная линия) и турбулентных (сплошная линия) фаз и От от интенсивности шума В. Величина В = 84 примерно соответствует случаю, когда амплитуда шума равна амплитуде исходного сигнала х(Ь)
4. Характер перемежающегося поведения в спонтанной неконвульсивной судорожной активности у крыс
Перейдем теперь от рассмотрения модельных систем к изучению поведения сложных реальных систем, демонстрирующих перемежающееся поведение. В качестве такого примера в данной работе приведены результаты анализа временной реализации записей электрической активности головного мозга крыс линии "Ав/Яу. Все крысы этой линии имеют генетическую предрасположенность к абсанс эпилепсии. Судорожная активность представлена пик-волновыми разрядами, которые возникают спонтанно, и не сопровождается конвульсиями [25].
В нашем эксперименте запись электрической активности мозга проводили с поверхности коры больших полушарий (электрокортикограмма, ЭКоГ) у интактных
свободно-подвижных крыс посредством хронически вживленных электродов. Частотный диапазон ЭКоГ составил 1-500 Гц, амплитуда сигнала варьировала от 25 до 800 мкВ. Более подробное описание метода приводится в работе [26].
Длительность ЭКоГ записей составляла от 6 часов до 4 суток. В работе были использованы как самцы (п = 5), так и самки (п = 6) вышеупомянутой линии крыс. Типичная ЭКоГ, регистрируемая в ходе наблюдений, представляла собой чередование низкоамплитудных (от 25 до 100 мкВ) полиритмических участков ЭКоГ, соответствующих «нормальному» функционированию головного мозга (будем называть
Рис. 10. a - Фрагмент электроэнцефалограммы, представляющий собой ламинарные фазы (нормальное поведение), перемежающиеся турбулентными всплесками (эпилептическая активность показана серыми прямоугольниками). б - Проекция вейвлет-поверхности \W(s, t)\, соответствующая электроэнцефалограмме. По оси абсцисс отложено время t, по оси ординат - временные масштабы s. Интенсивность градаций серого цвета пропорциональна абсолютной величине вейвлет-поверхности \ W(s, t)\. Сплошная линия ограничивает область влияния краевых эффектов, достоверные данные располагаются ниже этой линии. в - Зависимость энергии w(t), приходящейся на диапазон характерных временных масштабов (si, s2). Величина порога А = 0.5 показана пунктирной линией
такие участки ламинарными), и участков высокоамплитудной (от 100 до 800 мкВ) генерализованной активности с относительно стабильной несущей частотой от 8 до 12 Гц, соответствующих эпилептическим припадкам (условимся называть эти участки турбулентными фазами).
На рис. 10, а приведен фрагмент типичной ЭКоГ, состоящей из фоновой активности (ламинарных фаз), перемежающейся турбулентными всплесками, соответствующими эпилептическим пик-волновым разрядам (серые прямоугольники на рис. 10, а, в). Пик-волновые разряды были выделены экспертом по всей длине ЭКоГ на основании критериев, описанных в [25]. Проекция вейвлет-поверхности \Ш(в,Ь)\ показана на рис. 10, б. Отчетливо видно, что фрагменты вейвлет-поверхности, соответствующие ламинарным и турбулентным участкам, оказываются принципиально различными. Зависимость от времени энергии вейвлет-спектра, приходящейся на диапазон временных масштабов в £ (в1,в2) в = 70 мс, в2 = 110 мс) показана на рис. 10, в. При наступлении турбулентной фазы величина 'ш(Ь) превышает порог А = 0.5, что позволяет легко выделять различные типы поведения во временной реализации. Выбор диапазона временных масштабов (в1;в2) обусловлен структурой вейвлет-поверхности Ш(в, Ь) в моменты времени, соответствующие ламинарным и турбулентным участкам. Мгновенные распределения энергии вейвлет-спектра Е(в, Ь) \^0 для разных участков ЭКоГ показаны на рис. 11.
В результате проведенных исследований было установлено, что чередование судорожной и неэпилептической («нормального» функционирования мозга) активности может быть хорошо описано в рамках динамики оп-оГГперемежаемости [5,9,21]. Такое заключение сделано на основании анализа распределений длительностей фаз «нормального» функционирования мозга (длительностей ламинарных фаз), выделенных из ЭКоГ с помощью вышеописанной методики (см. также [6]).
На рис. 12 показаны распределения N(¿ь) длительностей ламинарных фаз для одной из крыс с учетом статистики по всем четырем дням, в течение которых проводились измерения. Распределения построены в двойном логарифмическом масштабе, точками нанесены экспериментальные данные по длительности ламинарных фаз,
Рис. 11. Типичные распределения энергии вейвлет-спектра по временным масштабам Е(в,10) для ламинарных (а) и турбулентных (б) фаз в различные моменты времени ¿0. Обратите внимание, что диапазоны значений величины Е(в) различаются на рис. (а) и (б) в десять раз
N 1
N1
10° -
101 -
102-
10°,
102п
ЮЦ
ю-1
' I 1111
ю-1
102
103
103
а
Рис. 12. Распределения N длительностей ламинарных фаз для одной из крыс (в двойном логарифмическом масштабе), точками (•) нанесены данные, полученные в ходе наблюдений. Прямые линии аппроксимируют экспериментально полученные распределения и соответствуют степенной зависимости с показателем степени а = -3/2: а - ночное время; б - дневное время
полученные в ходе наблюдений. Поскольку появление судорожной активности зависит от времени суток (наблюдается циркадная периодичность в частоте появления судорог [27]), анализ данных проводился отдельно для светлой и темной фазы светового режима. Рис. 12, а соответствует данным, полученным в ночное время, когда судорожная активность резко возрастает, поэтому статистика получается существенно более богатая). На рис. 12, б показано аналогичное распределение, построенное для судорожной активности, зарегистрированной в дневное время, когда число разрядов существенно меньше.
Из рис. 12 видно, что распределения длительностей ламинарных фаз, которые были получены в эксперименте (точки на рисунке), в двойном логарифмическом масштабе близки к прямой линии, характеризующей степенной закон распределения с показателем степени а = -3/2. Прямые сплошные линии на рис. 12, а, б аппроксимируют экспериментально полученные распределения и соответствуют степенной зависимости с показателем степени а = -3/2. Соответствующие параметры аппроксимации определялись численно путем минимизации среднеквадратичного отклонения теоретического распределения от экспериментальных данных.
Важно отметить, что перемежающееся поведение было диагностировано как для светлого, так и темного периодов, которые сильно различались по интенсивности возникновения судорог. При этом у данных животных в темную фазу число судорог оказывалось большим по сравнению со светлой фазой. Это привело, соответственно, к уменьшению средней длительности ламинарных фаз. В то же самое время, уменьшение среднего значения длительности ламинарных фаз подчинялось тому же степенному закону с показателем степени а = -3/2, следовательно, характер распределения длительностей ламинарных фаз не зависел от интенсивности турбулентных фаз.
Аналогичные результаты были получены при анализе судорожной активности данного типа на ЭКоГ всех экспериментальных животных (11 крыс). Таким образом, можно сделать вывод, что чередование судорожной и «нормальной» активности головного мозга крыс "Ав/Яу является перемежаемостью оп-оГГ типа и подчиняется универсальной закономерности, свойственной для данного типа перемежающегося поведения.
Заключение
В настоящей работе описан метод выделения ламинарных и турбулентных фаз во временных реализациях нелинейных систем, демонстрирующих перемежающееся поведение. Данный метод основан на непрерывном вейвлет-преобразовании и, как следует из данной работы, может быть с успехом применен для изучения различных типов перемежаемости в системах различной природы, включая живые системы (нейронные ансамбли, головной мозг [6,28], кардиоваскулярная система человека [29,30]). Предложенный метод является устойчивым к влиянию шумов и флуктуа-ций, искажающих исходную временную реализацию. Высокая эффективность предложенного метода обусловлена различающимися характерными временными масштабами, на которые приходится наибольшая доля энергии вейвлет-спектра, во время ламинарных и турбулентных фаз.
Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований (проекты 05-02-16273 и 07-02-00044), Программы поддержки ведущих научных школ РФ (НШ-4167.2006.2) и Программы Президиума РАН «Фундаментальные науки - медицине». Авторы благодарят также за финансовую поддержку Фонд некоммерческих программ «Династия».
Библиографический список
1. Берже П., Помо И., Видаль К. Порядок в хаосе. М.: Мир, 1991.
2. Manneville P., Pomeau Yv. Different ways to turbulence in dissipative dynamical systems // Physica D. 1980. Vol. 1(2). P. 167.
3. Rosenblum M.G., Pikovsky A.S., Kurths /.From phase to lag synchronization in coupled chaotic oscillators // Phys. Rev. Lett. 1997. Vol. 78(22). P. 4193.
4. Boccaletti S., Valladares D.L. Characterization of intermittent lag synchronization // Phys. Rev. E. 2000. Vol. 62(5). P. 7497.
5. Hramov A.E., Koronovskii A.A. Intermittent generalized synchronization in unidirectionally coupled chaotic oscillators // Europhysics Letters. 2005. Vol. 70(2). P. 169.
6. Короновский А.А., Кузнецова Г.Д., Мидзяновская И.С., Ситникова Е.Ю., Тру-бецков Д.И., Храмов А.Е. Закономерности перемежающегося поведения в спонтанной неконвульсивной судорожной активности у крыс // Доклады Академии Наук, 2006.
7. Berge P., Pomeau Y., Vidal Ch. L'Ordre Dans Le Chaos, 1988.
8. Dubois M., Rubio M., Berge P. Experimental evidence of intermittencies associated with a subharmonic bifurcation // Phys. Rev. Lett. 1983. Vol. 51. P. 1446.
9. PlattN., Spiegel E.A., Tresser C. On-off intermittency: a mechanism for bursting // Phys. Rev. Lett. 1993. Vol. 70(3). P. 279.
10. Pikovsky A., Osipov G., Rosenblum M., Zaks M., Kurths /.Attractor-repeller collision and eyelet intermittency at the transition to phase synchronization // Phys. Rev. Lett. 1997. 79(1). P. 47.
11. Zhan M, Wei G.W., Lai C.H. Transition from intermittency to periodicity in lag synchronization in coupled Rossler oscillators // Phys. Rev. E. 2002. Vol. 65(3). 036202.
12. Короновский А.А., Храмов А.Е. Об эффективном анализе перехода к хаосу через перемежаемость с помощью вейвлетного преобразования // Письма в ЖТФ. 2001. Т. 27(1). С. 3.
13. Короновский А.А., Храмов А.Е. Непрерывный вейвлетный анализ и его приложения // М.: Физматлит, 2003.
14. Daubechies I. Ten lectures on wavelets // SIAM, 1992.
15. Torresani B. Continuous wavelet transform//Paris: Savoire, 1995.
16. Короновский А.А., Тыщенко А.А., Храмов А.Е. Исследование распределения турбулентных фаз при разрушении синхронизации с запаздыванием // Письма в ЖТФ. 2005. Т. 31(21).
17. Hramov A.E., Koronovskii A.A. An approach to chaotic synchronization // Chaos. 2004. Vol. 14(3). P. 603.
18. Шустер Г. Детерминированный хаос. М.: Мир, 1988.
19. Pecora L.M., Carroll T.L. Driving systems with chaotic signals // Phys. Rev. A. 1991. Vol. 44(4). P. 2374.
20. Murali K., Lakshmanan M.Drive-response scenario of chaos syncronization in identical nonlinear systems // Phys. Rev. E. 1994. Vol. 49(6). P. 4882.
21. Heagy J.F., Platt N., Hammel S.M. Characterization of on-off intermittency // Phys. Rev. E. 1994. Vol. 49(2). P. 1140.
22. Короновский А.А., Москаленко О.И., Храмов А.Е. Новый тип универсальности при хаотической синхронизации динамических систем // Письма в ЖЭТФ. 2004. Т. 80(1). С. 25.
23. Hramov A.E., Koronovskii A.A., Kurovskaya M.K., Moskalenko O.I. Synchronization of spectral components and its regularities in chaotic dynamical systems. Phys. Rev. E. 2005. Vol. 71(5). P 056204.
24. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.T. Numerical Recipes // Cambridge University Press, Cambridge, 1997.
25. Van Luijtelaar E.L., Coenen A.M. Two types of electrocortical paroxysms in an inbred strain of rats // Neurosci Lett. 1986. Vol. 70(3). P. 393.
26. Мидзяновская И.С. Два типа разрядов «пик-волна» на электроэнцефалограмме крыс линии wag/rij, генетической модели absence эпилепсии // Журн. высш. нерв. деят. 1999. Т. 49(5). С. 855.
27. Luijtelaar E.L, Coenen A.M. Circadian rhythmicity in absence epilepsy in rats // Epilepsy Res. 1998. Vol. 2(5). P. 331.
28. Tass P.A. et al. Synchronization tomography: A method for three-dimensional localization of phase synchronized neuronal populations in the human brain using magnetoencephalography. Phys. Rev. Lett. 2003. Vol. 90(8). P 088101.
29. Hramov A.E., Koronovskii A.A., Ponomarenko V.I., Prokhorov M.D. Detecting synchronization of self-sustained oscillators by external driving with varying frequency // Phys. Rev. E. 2006. Vol. 73. P. 026208.
30. Короновский А.А., Пономаренко В.И., Прохоров М.Д., Храмов А.Е. Изучение синхронизации автоколебаний по унивариантным данным при изменении частоты внешнего воздействия с использованием вейвлетного анализа // Письма в ЖТФ. 2006. Т. 32(11). С. 81.
Саратовский государственный Поступила в редакцию 20.03.2007
университет После доработки 7.06.2007
Институт высшей нервной деятельности и нейрофизиологии РАН
NICI-Biological Psychology, Radboud University Nijmegen, Nijmegen, the Netherlands
APPLICATION OF CONTINUOUS WAVELET TRANSFORM TO ANALYSIS OF INTERMITTENT BEHAVIOR
A.A. Koronovskii, I.M. Minyukhin, A.A. Tyshchenko, A.E. Hramov, I.S. Midzyanovskaya, E.Yu. Sitnikova, G. van Luijtelaar, C.M. van Rijn
Effective method of signals analysis based on the continuous wavelet transform is proposed in this paper. Application of this method to estimation of mean value both of laminar and turbulent phase durations corresponding to different types of intermittent behavior is considered including analysis of time series produced by living systems. It is shown that the proposed method is stable to noise and fluctuations distorting the initial time series.
Минюхин Игорь Михайлович - родился 24 января 1987 года. Студент 3 курса факультета нелинейных процессов Саратовского государственного университета. Область научных интересов: применение методов нелинейной динамики к изучению систем различной природы, в том числе социальных, биологических, геологических и т.п. Участник студенческих конференций, автор публикации в сборнике докладов «Нелинейные дни-2005».
Тыщенко Александр Александрович - родился в Саратове (1984), окончил Саратовский государственный университет, факультет нелинейных процессов (2006), продолжает обучение в магистратуре по направлению радиофизика. Область научных интересов - исследование сложных режимов нелинейных динамических систем. Имеет одну публикацию по данному направлению.
Мидзяновская Инна Станиславовна - родилась в Московской области (1974), окончила физический факультет (кафедра биофизики) Московского государственного университета им М.В. Ломоносова в 1997 году. После окончания МГУ работает в ИВНД и НФ РАН. Победитель конкурса Молодых Ученых на 23-м международном эпилептологическом конгрессе в Праге, Чешская Республика (1999) и на 1-м симпозиуме по нейроповеденческой фармакологии в Окаяме, Япония (2004), финалист конкурса Молодых Ученых на 32-й ежегодной встрече Европейского Гистаминового Сообщества, Ноордвикерха-ут, Нидерланды (2000). Защитила диссертацию «Абсансная и смешанная формы эпилепсии у крыс WAG/Rij: характеристики и аминергическая модуляция» в Университете Радбауд г. Наймегена (Нидерланды, 2006). Автор 20 статей по экспериментальной эпилептологии, нейрохимии, анализу ЭЭГ. Рецензент журналов «Neuroscience Letters», «Pharmacology, Biochemistry and Behavior». E-mail:[email protected]
Ситникова Евгения Юрьевна - родилась в Волгограде (1974). Окончила ветеринарно-биологический факультет Московской академии ветеринарной медицины и биотехнологии по специальности «биофизика» (1996), кандидат биологических наук, научный сотрудник лаборатории нейроонтогенеза Института высшей нервной деятельности и нейрофизиологии РАН, Москва. Научные интересы в области нейрофизиологии включают математические методы анализа активности головного мозга человека и животных, изучение влияния факторов окружающей среды на формирование функций мозга, адаптивные механизмы поведения и памяти; морфо-функциональные механизмы взаимодействий нервных клеток в ходе индивидуального развития в норме и патологии. Автор более 20 научных трудов в отечественной и зарубежной печати.
G. van Luijtelaar is the senior researcher at the NICI, Dept of Biological Psychology. He received his PhD from Radboud University Nijmegen. He studied Biological Psychology in Tilburg University. His research interests are currently «mechanisms of normal pathological thalamo-cortical interactions». E-mail:[email protected]
C.M. van Rijn - received her bachelors of chemistry and in 1984 she graduated as a medical doctor from the University of Leiden, the Netherlands. In 1989 she received her PhD in Nijmegen, the Netherlands on complementary in vivo (EEG) and in vitro studies (receptor binding) concerning the interrelationship of GABA antagonists and epilepsy. Since 1988 she is a docent in NICI (Nijmegen Institute for Cognition and Information), Department of Biological Psychology, Radboud University of Nijmegen. Her research is concentrated on the influence of GABAA -ergic drugs on the EEG of rats and of humans. Clinical topics of her special interest are epileptic phenomena, anaesthesia and neuropathic pain. E-mail:[email protected]
