CC BY
21
Нелинейная
динамика и неиронаука
Изв. вузов «ПНД», т. 21, № 5, 2013 УДК 621.391.01
ПОСЛЕДОВАТЕЛЬНАЯ ПЕРЕКЛЮЧАТЕЛЬНАЯ АКТИВНОСТЬ В АНСАМБЛЕ НЕИДЕНТИЧНЫХ СИСТЕМ ПУАНКАРЕ
А. О. Михайлов, М. А. Комаров, Г. В. Осипов
В работе анализируется последовательная переключательная активность в ансамбле тормозно связанных систем Пуанкаре. В ходе исследования обнаружено существование гетероклинического контура в фазовом пространстве при асимметричной тормозной связи при определённом значении параметра скорости связи.
Изучена динамика ансамбля тормозно и диффузионно связанных неидентичных систем Пуанкаре. Построены приближенные бифуркационные диаграммы для всех качественно различных режимов активности сети и определены области пространства параметров, соответствующие различным динамическим режимам, таким как мультистабиль-ность, вымирание, модуляция, последовательная пачечная активность и синхронизация.
Ключевые слова: Последовательная активность, гетероклинический контур, синхронизация.
Введение
Многие динамические процессы в нейронных ансамблях представляют собой последовательные переключения активности между отдельными элементами и (или) группами элементов. Такая последовательная нейронная активность может быть связана со многими физиологическими функциями нервной системы [1-3]. Генерация и распространение возбуждения между отдельными элементами, а также между группами нейронов играют важнейшую роль в функционировании мозга и нервной системы в целом. Например, известно, что подобный тип активности является неотъемлемым свойством нейронных сетей в сенсорных и двигательных системах животных [4]. Также, определённый отдел мозга птиц генерирует последовательности пачечной активности, которые, в свою очередь, управляют голосовым аппаратом и обеспечивают воспроизведение песни [5].
С точки зрения нелинейной динамики причиной последовательной активности в нейронных сетях может быть существование устойчивого гетероклинического контура в фазовом пространстве динамической системы, моделирующей активность
нейронной сети [6,7]. В основе генерации последовательных переключений активности лежит принцип беспобедительной конкуренции [8]. Суть этого принципа заключается в последовательном прохождении изображающей точки в окрестности определённой седловой траектории, что означает активацию определённого нейрона или группы нейронов [9,10]. Таким образом, устойчивый гетероклинический контур в фазовом пространстве является математическим образом последовательных переключений активности в ансамблях связанных нейронов.
В данной работе исследуется вопрос о существовании устойчивого гетерокли-нического контура между седловыми предельными циклами в ансамбле связанных систем Пуанкаре. Данная система была выбрана в связи с тем, что она допускает аналитическое исследование, в отличие от моделей типа Ходжкина-Хаксли, сложность которых приводит к тому, что изучение их динамики носит преимущественно численный характер. Описанные далее эффекты, обнаруженные для рассматриваемой системы, аналогичны полученным ранее для модели Ходжкина-Хаксли [11], что подтверждает их универсальный характер.
1. Модель
Рассмотрим следующую модель тормозно связанных систем Пуанкаре
Xг = -Уг + Хг[А(вг) - (х2 + у2)],
уг = хг + Уг[А(вг) - (х2 + у2)\, N
Твг = ^ ЯЦ Е (Р2) 3=1,3=
- вг
г =
(1)
где
Рг = ^х1 + Уг2,
Т » 1,
г = 1,Ы,
А(вг) = 1 - в2, г = 1,^, Е (Р2) =_1___1_
1+ е-(р1-хо)/к 1 + ехо/к>
г = 1,ы.
Определим возможный физический смысл переменных, входящих в (1): хг - мембранный потенциал г-й клетки; уг - активационная переменная, отвечающая за совокупность ионных токов; вг - переменная, отвечающая за тормозные связи г-го элемента; т - параметр скорости связи; дгз - параметр, отражающий силу влияния тормозной связи ]-й клетки на г-ю; N - число элементов в сети; Е - активационная функция клетки; хо и к - параметры активационной функции Е. В фазовом пространстве изолированного элемента системы (1) существует устойчивый предельный цикл, который представляет собой окружность радиуса 1. Переходя к полярной
системе координат, можно переписать (1) в виде
' рг = рг[А(вг) - р2],
Ф г = 1,
N
Твг
Е За Р (р.2)
- вг
з=1, з=г
{г = 1,Ж.
Рассматривая амплитуды, можно получить, что предельным циклам в системе координат (хг, у г, в г) соответствуют состояния равновесия в системе координат (рг, вг).
Для того, чтобы понять динамику больших (Ж > 3) систем, исследуем сначала динамику двух связанных систем Пуанкаре, рассматривая двумерное пространство амплитуд (р!, р2). Полагая р^ =0, ] = 1, 2, получаем
' (3! = р![А(в!) - р.], Тв?! = З!2Р(р2) - в!, р2 = р2[А(в2) - р2], Тв2 = З2!Р(р?) - в2.
(3)
Область медленных движений в (3) имеет вид
( р?[А(в!) - р?] = 0, 1 р2[А(в2) - р2] =0.
Рассмотрим только одну подобласть области медленных движений
р!
р22
1 - в2,, в2 < 1, в? > 1,
1 - в2, в2 < 1,
в2 > 1.
(4)
В области быстрых движений имеем
в? ^ еош^ в2 ^ еош!.
После подстановки (4) в (3) получаем
Тв1
Тв2
Я12
Я12
1 + е-{1-822-х0)/к 1 + ехо/к
д21
д21
1+ е-(1-«1-х0)/к 1+ ехо/к
в1 < 1,
кв2 < 1.
в1,
в2,
2. Симметричная связь
Рассмотрим случай симметричных больших связей, то есть положим
Я12 = Я21 = д> 2. (6)
Покажем, что в (5) при (6) нет замкнутых контуров, состояний равновесия и что все траектории, находясь в области {в1;в2 € (0; 1) & в1 = в2}, выходят из нее. Будем рассматривать не весь квадрат (0; 1) х (0; 1), а лишь его половину, то есть треугольник{
<в1 < 1 (7)
\в1 <в2 < 1.
Это правомерно, так как система (5) инвариантна относительно замены в1 на вг и наоборот. Пусть к ^ 1. Полагая в1 = 0, из (5) получаем, что в1 > 0, а значит, траектории выходят из отрезка
{ в1 = 0, | 0 < вг < 1.
Если в1 = вг, нетрудно прийти к выводу о том, что траектории из треугольника (5) не могут пересекать отрезок {
I в1 = в2, | 0 < в2 < 1.
Применяя критерий Бендиксона, получаем, что замкнутых контуров в треугольнике (7) не существует.
Покажем, что в (7) нет состояний равновесия. Уравнения для отыскания состояний равновесия в (5) имеют вид
в1 =
д
д
1 + е-(1-822-хо)/к 1 + ехо/к'
(8)
в2
1 + е-(1-8^-х о)/к 1 + ехо/к'
Подставляя во второе уравнение системы (8) значение в^, выраженное из первого
д
д
уравнения той же системы, получаем, что состояний равновесия системы (5) в области {«!, в2 € (0; 1) & = в2} нет.
Таким образом, показано, что динамика (3) в области (4) такова, что все траектории, находящиеся в квадрате [0; 1] х [0; 1], остаются там. Тогда, в силу (4), получаем, что изображающая точка системы (3) приходит либо в точку (0,1), либо в точку (1, 0). Данный факт проиллюстрирован на рис. 1.
Р2 1.0
0.8
0.6
0.4
0.2
0
N
1 ж
¿¿о/ ->-II-L
0 0.2 0.4 0.6 0.8 1.0 р Рис. 1. Фазовый портрет системы (3) в двумерном пространстве амплитуд. Параметры: т = 100, к = 0.01, хо = 0.25, 012 = 921 =3
3. Асимметричная связь
Рассмотрим систему (3) при асимметричных связях, например, при следующих значениях параметров:
gi2 > 2,
521 < 1
(9)
Тогда #21^(р1) € [0; #21], что, в свою очередь, в силу (3) означает, что в2 € € [0; д21]. А значит, в силу (3), Р2 € [1 - д^; 1] и F(р2) > 1/2, откуда следует, что g12F(р2) > 1. Это заключение позволяет сделать вывод о том, что при £ ^ «1 > 1, поэтому А(«1) = = 1 — в2 < 0. В итоге, при £ ^
Г Р1 ^ 0, \ Р2 ^ 1
Таким образом, мы показали существование гетероклинического контура при асимметричных связях. В плоскости (Р1, Р2) динамика системы (3) представлена на рис. 2.
Рассматривая отдельно динамику двух связанных систем Пуанкаре в плоскостях (рг, pj), где г = ] и г,] = 1,3, а
Рис. 2. Фазовый портрет системы (3) в плоскости (pi, p2). Параметры: т = 100, k = 0.01, x0 = 0.25,
gi2 = 3, g2i = 0.5
Рис. 3. Фазовый портрет (2) при N = 3. Параметры: т = 100, k = 0.01, xo = 0.25, gi2 = gsi =
= g23 = 4, gi3 = g2i = g32 = 0.5
после этого комбинируя полученные результаты в пространстве (р1, р2, Рз), получаем динамику (2) (при N = 3) в целом, которая изображена на рис. 3. На данном рисунке изображены три гетероклинические орбиты между состояниями (1) О1 и О2, (и) 02 и 03, (Ш) 03 и 01.
4. Бифуркационный анализ
Рассмотрим вопрос о бифуркации, которая происходит при изменении параметров тормозных связей д^. Для этого отыщем состояния равновесия системы (3). Они определяются из уравнений
(Р1(1 - д!гр2(р2) - р1) = 0, рг(1 - д\х?2(р2) - р2) = 0,
в1 = д12^ (Р2),
в2 = д21^ (Р2)-Рассмотрим состояние равновесия системы (10)
л л „2 т^2/л\
(10)
Р1 = 1 - д12^ (р2),
Л л Л Т^2/Л\
(11)
р2 = 1 - д21^ (р1).
Систему (11) можно представить в следующей форме:
р1 = 1 - д2^ 2[1 - д22^ 2(р1)], (12)
Корни уравнения (12) являются нулями следующей функции:
н(Р1) = 1 -д212р2[1 - 2(р1 )] - р1.
Будем искать состояние равновесия (12) численно. Положим, что т = 100; к = 0.01; х0 = 0.25; д12 = 3. Значение д21 будем изменять от ц = 0.88 < 1 до нуля. Функция Н (р1) в зависимости от значения д21 имеет два, один или нуль корней (см. фрагменты а, в, д на рис. 4). Динамика в плоскости (х2,у2) в этом случае показана на фрагментах б, г и д на рис. 4. Таким образом, имеем седло-узловую бифуркацию предельного цикла.
Динамика (3) при выбранных значениях параметров и д21 = 0.88 показана на рис. 5.
Рассмотрим вопрос устойчивости гетероклинического контура. Для этого исследуем (2) при N = 3. Принимая во внимание условие (9), вычислим значения седловых величин, относящихся к входящим в гетероклинический контур состояниям равновесия. Положим, д1з = дз2 = д21 = д. Используя теорему об устойчивости контура между седловыми предельными циклами из [9], получим условие устойчивости гетероклинического контура в зависимости от отношения двух параметров -т и д2.
1
(1 - д2Е2(1))3т3'
Гетероклинический контур устойчив тогда и только тогда, когда V > 1. Найдём условие устойчивости. Пусть V =1. Если к = 0.01, хо = 0.25, то Е2 (1) = 1. Тогда имеем
1
1 - д2
т
Таким образом, получаем кривую, которая выделяет область устойчивости ге-тероклинического контура и область существования предельного цикла в пространстве (р1, р2, р3) в окрестности контура. Режим переключательной активности переменной х рассматриваемой нами сети при N = 3 изображён на рис. 6.
Рис. 4. Динамика системы (1) в плоскости (х2 , У2). Имеем один устойчивый и один неустойчивый предельный цикл внутри окружности радиуса 1 плоскости (х2 ,У2) в случае (а-б), полуустойчивый предельный цикл в случае (в-г) и не имеем предельных циклов в случае (д-е). Параметры:
а-б - 921 = 0.88, в-г - 921 = 0.8718206, д-е - 921 = 0.8
Рис. 5. Фазовый портрет (3) при т = 100, к = Рис. 6. Граница области устойчивости гетерокли-= 0.01, хо = 0.25, 912 = 3 и 921 = 0.88 нического контура и динамика переменной х
5. Динамика ансамбля диффузионно связанных неидентичных систем Пуанкаре
Известно, что помимо тормозной связи клетки, составляющие сеть, связаны между собой и электрически [4]. Электрическая связь между нейронными клетками моделируется в дифференциальных уравнениях через центральную разность второго порядка, которая физически работает как диффузионная связь.
Рассмотрим модель связанных систем Пуанкаре, обладающих тормозной и диффузионной связью
N
Xг = -WiУi + Хг[Л(вг) - (X2 + у2)] + й Е X - Хг),
з=1, з=г
N
Уг = 'ШгХг + уг[А(вг) - (х2 + у2)] + й ^ (Уз - Уг)>
< з=1, з=г
N
хёг = Е дгзР(Р2) -
з=1, з=г
г = ГЖ
где параметр й можно интерпретировать как силу электрической связи между элементами сети. Изменяя частотную расстройку между элементами при N = 3, будем изучать влияние параметра диффузионной связи й. Частотную расстройку будем выбирать особым образом, для чего введём параметр А = W2 - wl = w3 - W2, который будет отражать степень неидентичности осцилляторов.
Рассмотрим случай больших симметричных связей. Приближенная бифуркационная диаграмма для рассматриваемого ниже случая и иллюстрации к существующим при этом режимам приведены на рис. 7. При малых значениях А (приблизительно до 0.6) какой-то один из элементов сети (в зависимости от начальных условий) со временем подавляет активность двух других, то есть наблюдается мультистабиль-ность. В то же время, было установлено, что параметр диффузионной связи оказывает сильное влияние на частотную динамику элементов сети в данной области. Так, область 4 характеризуется одинаковой для всех элементов частотой, равной частоте подавляющего элемента. В то время как в области 5, частоты всех элементов с увеличением диффузионной связи стремятся к средней частоте. Частотный переход из области 4 в область 5 является мягким. Дальнейшее увеличение частотной расстройки приводит к возникновению синхронизации наряду с мультистабильностью. Частотная динамика области 2 такова, что все элементы сети синхронизируются на средней частоте. Заметим, что частотный переход из области 4 в область 2 является жёстким. При больших значениях А (приблизительно с 1.3), наряду с указанными режимами, наблюдается ещё режим модуляции, который характеризуется достаточно большим разбросом по частоте и амплитуде колебаний переменной Х. Дальнейшее
увеличение неидентичности элементов сети (примерно с А = 1.78) приводит к появлению ещё одного режима активности сети - к так называемому «вымиранию», то есть режиму, при котором все элементы находятся в состоянии покоя, никак не проявляя свою активность.
Рассмотрим случай асимметричных связей, при котором существует предельный цикл в системе (1) в окрестности гетероклинического контура. Рассматривая динамику сети при малых значениях А (примерно до 0.6), имеем последовательную пачечную активность с постоянными для всех элементов длинами пачек. При изменении параметра диффузионной связи частота пачек и их профиль претерпева-
Рис. 7. а - бифуркационная диаграмма для параметров: т = 100, к = 0.01, х0 = 0.25, 912 = 931 = 923 = 913 = 921 = 932 = 4, А = ш2 — ш1 = ш3 — ш2. Режимы: 1 - вымирание, 2 - синхронизация (случай б), 3 - модуляция (случай в), 4-5 - мультистабильность (случаи г-д)
ют изменения. Увеличивая силу диффузионной связи, можно заметить уменьшение частоты пачечной активности. Это связано с увеличением времени готовности элемента стать активным, при возрастании й. Фактически наблюдается увеличение порога активности элемента при возрастании силы диффузионной связи. Дальнейшее увеличение частотной расстройки (примерно с А = 0.6) приводит к возникновению синхронизации наряду с последовательной пачечной активностью . Однако, начиная приблизительно с А = 0.7, возникает модуляция, которая является переходным режимом от последовательной пачечной активности к синхронизации. При больших значениях уровня неидентичности элементов сети (приблизительно с А = 1.78), наряду с указанными режимами наблюдается режим вымирания. Приближенная бифуркационная диаграмма в этом случае и иллюстрации к существующим при этом режимам показаны на рис. 8.
Рассмотрим случай асимметричных связей, при котором существует гетеро-клинический контур в системе (1). При малых значениях А (примерно до 0.45), имеем режим модуляции при всех интересующих нас значениях диффузионной связи. Продолжая увеличение частотной расстройки, получаем, что наряду с модуляцией
Рис. 8. а - бифуркационная диаграмма в случае следующих значений параметров: т = 100; к = 0.01; хо = 0.25; §12 = 931 = 923 = 4; 913 = д21 = 932 = 0.5; А = т — т = тз — т2. Режимы: 1 - вымирание, 2 - синхронизация (случай б), 3 - модуляция (случай в), 4 - последовательная пачечная активность (случай г)
д 1.6
1.2
0.8
0.4
а о 0.2 0.4 0.6 0.8 й
1.0 0.5 X! О 0.5 -1.0 1.0 0.5 х2 0 0.5 -1.0 1.0 0.5 х3 0 0.5
-1.0
О
б
Рис. 9. а - бифуркационная диаграмма в случае следующих значений параметров: т = 0.01; к = 0.01; Хо = 0.25; §!2 = дз1 = 923 = 4; 913 = д21 = дз2 = 0.5; А = -Ш2 - "Ш1 = "Шз - Ы2- Режимы: 1 - вымирание, 2 - синхронизация (случай в), 3 - модуляция (случай б)
существует синхронизация. При больших значениях А (начиная, приблизительно с 1.78), одновременно с указанными режимами наблюдается режим вымирания. Приближенную бифуркационную диаграмму для этого случая и иллюстрации к существующим при этом режимам можно увидеть на рис. 9.
Заключение
В статье проведено исследование ансамбля тормозно связанных систем Пуанкаре. Подобную сеть феноменологически можно рассматривать как нейродинамиче-скую модель. Тем не менее, она демонстрирует основные эффекты, наблюдающиеся в биологически реалистичных моделях типа Ходжкина-Хаксли. Было показано, что при определённых значениях параметров существует устойчивый гетероклини-ческий контур. Так как гетероклинический контур является математическим образом последовательной активности в нейронных ансамблях, то, по сути, был выявлен данный режим активности у определённого класса нейронов. Таким образом, используя полученные результаты, возможно воссоздание режима последовательной активности в реальных биологических нейросетях и дальнейшее исследование этого эффекта, чья роль в нервной системе обсуждалась выше.
Также было проведено численное исследование динамики переменной х, которую можно интерпретировать как мембранный потенциал, с частотной расстрой-
кой осциллирующих элементов, учитывая, помимо тормозной связи, диффузионную. В результате данного исследования построены приближенные бифуркационные диаграммы для всех качественно различных случаев динамики сети. Определены области пространства параметров, соответствующие различным динамическим режимам, таким как мультистабильность, вымирание, модуляция, последовательная пачечная активность и синхронизация.
Работа выполнена при финансовой поддержке РФФИ (гранты 10-02-00940, 11-02-92003, 11-07-97013) и ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 гг. (№ контракта 02.740.11.0919).
Библиографический список
1. KomarovM.A., Osipov G.V., Suykens J.A.K. andRabinovich M.I. Numerical studies of slow rhythms emergence in neural microcircuits: Bifurcations and stability // CHAOS. 2009. Vol. 19. 015107.
2. Galan R., Sasche S., Galicia C.G. and Herz A.V. Odor-driven attractor dynamics in the antennal lobe allow for simple and rapid olfactory pattern classification // Neur. Comput. 2004. Vol. 16. P. 999.
3. LeviR., VaronaP., Arshavsky Y.I., Rabinovich M.I. and Selverstone A.I. Dual sensory-motor function for a molluskan statocyst network // J. Neurophysiol. 2004. Vol. 91. P. 336.
4. Rabinovich M.I., Varona P., Selverston A.I. and Abarbanel H.D.I. Dynamical principles in neuroscience // Rev. Mod. Phys. 2006. Vol. 19. 015107.
5. Hahnloser R.H.R., Kozhevnikov A.A. and Fee M.S. An ultra-sparse code underlies the generation of neural sequences in a songbird // Nature (London). 2002. Vol. 419. P. 65.
6. Ashwin P., Burylko O. and Maistrenko Y Bifurcation to heteroclinic cycles and sensitivity in three and four coupled phase oscillators // Physica D. 2008. Vol. 237. P. 454.
7. Ashwin P. and Field M. Heteroclinic networks in coupled cell systems // Archive for Rational Mechanics and Analysis. 1999. Vol. 148. P. 107.
8. Seliger P., Tsimring L.S. and Rabinovich M.I. Dynamics-based sequential memory: Winnerless competition of pattern // Phys. Rev. E. 2003. Vol. 67. 011905.
9. Afraimovich VS., Rabinovich M.I. and Varona P. Heteroclinic contours in neural ensembles and the winnerless competition principle // Int. J. Bif. and Chaos. 2004. Vol. 14. 1195.
10. Afraimovich VS., Zhigulin V.P. and Rabinovich M.I. On the origin of reproducible sequential activity in neural circuits // CHAOS. 2004. Vol. 14. 1123.
11. Komarov M.A., Osipov G.V. and Suykens J.A.K. Sequentially activated groups in neural networks // Europhys. Lett. 2009. Vol. 86. 60006.
Нижегородский государственный Поступила в редакцию 25.04.2013
университет им. Н.И. Лобачевского После доработки 8.07.2013
SEQUENTIAL SWITCHING ACTIVITY IN THE ENSEMBLE OF NONIDENTICAL POINCARE SYSTEMS
A. O. Mikhaylov, M.A. Komarov, G. V. Osipov
Switching activity in the ensemble of inhibitory coupled Poicare systems is considered. The existence of heteroclinic contour in the phase space at the certain domain of parameter space has shown.
Dynamics of the ensemble of non-identical inhibitory and diffusively coupled systems of Poincare is considered. The approximate bifurcation diagrams for all qualitatively different regimes of the network activity have shown. There are areas of the parameter space corresponding to different dynamic regimes, such as multistability, extinction, modulation, bursting and synchronization.
Keywords: Sequential activity, heteroclinic contour, synchronization.
Михайлов Алексей Олегович - родился в 1989 году, окончил Нижегородский государственный университет им. Н.И. Лобачевского (2012). Аспирант кафедры теории управления и динамики машин факультета ВМК ННГУ. Область научных интересов: нейродинамика, синхронизация и последовательная активность в осцилляторных сетях, циклический хаос. Имеет 2 публикации (в соавторстве), в том числе 1 журнальную.
603950 Россия, Нижний Новгород, пр. Гагарина, 23 Нижегородский государственный университет им. Н.И. Лобачевского E-mail: aomikhailov@gmail.com
Комаров Максим Андреевич - родился в 1985 году, окончил Нижегородский государственный университет им. Н.И. Лобачевского (2008). В 2011 году защитил диссертацию на соискание ученой степени кандидата физико-математических наук по радиофизике. С 2013 года работает в университете Потсдама. Имеет 15 публикаций (в соавторстве).
603950 Россия, Нижний Новгород, пр. Гагарина, 23 Нижегородский государственный университет им. Н.И. Лобачевского E-mail: maxim.a.komarov@gmail.com
Осипов Григорий Владимирович - родился в Нижнем Новгороде (1960), окончил Нижегородский государственный университет им. Н.И. Лобачевского (1982). Защитил диссертацию по математическому моделированию на соискание ученой степени кандидата физико-математических наук в Научном совете по комплексной проблеме «Кибернетика» (1988) и доктора физико-математических наук по радиофизике (2004, ННГУ). С 1988 года работает в ННГУ, с 2007 года в качестве заведующего кафедрой «Теории управления и динамики машин». Соавтор монографий «Устойчивость. Структуры и хаос в нелинейных сетях синхронизации» и «Synchronization in oscillatory networks». Опубликовал более 120 научных статей (в том числе 2 обзора) по теории колебаний и волн и математическому моделированию.
603950 Россия, Нижний Новгород, пр. Гагарина, 23 Нижегородский государственный университет им. Н.И. Лобачевского E-mail: grosipov@gmail.com
