CC BY
21
ПОСЛЕДОВАТЕЛЬНАЯ АКТИВНОСТЬ В НЕЙРОННОМ АНСАМБЛЕ С ВОЗБУЖДАЮЩИМИ СВЯЗЯМИ
А. Г. Коротков, Г. В. Осипов
В статье предложена новая модель ансамбля нейроноподобных элементов, для моделирования которой используется обобщённая модель Лотки-Вольтерры с возбуждающими связями. Работа мотивирована тем, что возбуждающие связи составляют преобладающий тип взаимодействия между нейронами головного мозга. В работе показано, что в таком ансамбле в зависимости от связей между элементами существуют 2 режима: режим с устойчивым гетероклиническим циклом и режим с устойчивым предельным циклом.
Ключевые слова: Нейронный ансамбль с возбуждающими связями, модель Лотки-Воль-терры, численное моделирование.
Введение
Многие динамические процессы в нейронных ансамблях представляют собой последовательные переключения активности между отдельными элементами и (или) группами элементов. Такая последовательная нейронная активность может быть связана со многими физиологическими функциями нервной системы. Генерация и распространение последовательностей возбуждения между отдельными элементами, а также между группами нейронов играют важнейшую роль в функционировании мозга и нервной системы в целом. Например, известно, что подобный тип активности является неотъемлемым свойством нейронных сетей в сенсорных и двигательных системах животных [1]. Также, определенный отдел мозга птиц генерирует последовательности пачечной активности, которые, в свою очередь, управляют голосовым аппаратом и обеспечивают воспроизведение песни [2]. Существует ещё множество примеров, где ключевую роль играют последовательности возбуждения в нейронных сетях [3,4].
Чрезвычайно актуально рассмотреть последовательную активность в нейронных сетях с точки зрения нелинейной динамики. Существует гипотеза [1], что причиной последовательной активности в нейронных сетях может быть образование
устойчивого гетероклинического контура в фазовом пространстве динамической системы, моделирующей активность нейронной сети [5-8]. В основе генерации последовательных переключений активности лежит принцип winnerless competition (далее WLC), который можно перевести как принцип беспобедительной конкуренции [9]. Суть этого принципа как раз и заключается в существовании в фазовом пространстве устойчивого гетероклинического контура между особыми траекториями седло-вого типа (седловыми состояниями равновесия, седловыми предельными циклами и т.д.) [10,11]. Прохождение изображающей точки в окрестности определенной седло-вой траектории означает активацию определенного нейрона или группы нейронов. Таким образом, устойчивый гетероклинический контур в фазовом пространстве является математическим образом последовательных переключений активности в ансамблях связанных нейронов. Необходимым условием существования такого поведения является наличие ингибиторных связей между нейронами.
Динамические системы, моделирующие нейронную активность, делятся на две группы: биологически релевантные модели (модели на основе модели Ходжкина-Хаксли) и феноменологические (системы ФитцХью-Нагумо, Хиндмарш-Розе, ван дер Поля, накопление-сброс, Вилсона-Коуэна, Лотки-Вольтерры и др.). В работе [12] для моделирования нейронного ансамбля была использована обобщённая модель Лотки-Вольтерры. В работе [10] было показано, что в ансамбле ингибиторно связанных элементов Лотки-Вольтерры существует устойчивый гетероклинический цикл. В данной работе мы исследуем вопрос о переключательной активности в ансамбле возбуждающе связанных элементов Лотки-Вольтерры.
1. Модель ансамбля
Рассмотрим ансамбль из N элементов типа Лотки-Вольтерры с возбуждающими связями. Каждый элемент задаётся уравнениями
Рз* = (3 + ¿3)Рз(-1 + 3/Мз - Рз2), (1)
где р3- - уровень спайковой активности ]-го нейрона [1],
N
3 =Е 9г3Р(Рг) (2)
г=1
- входное воздействие на ]-й элемент со стороны остальных элементов; д^ - коэффициенты связей между г-м и ] -м элементами;
N
¿3 =Е д^Р(Рг) (3)
г=1
- дополнительный параметр, который равен 2 тогда и только тогда, когда активен элемент, который будет активироваться после j-го элемента;
( 1, Р > Я,
Р (Р) = { (4)
[0, р <д
- функция Хевисайда с порогом д = 0.999. Условимся считать элемент г активным, если Р(рг) = 1;
м ( Е Р (рг), 13 =0,
Ш3 = < г=1,дц =0 (5)
I 1 ^ =0
- число активных элементов, оказывающих воздействие на ] -й элемент, если таковых нет, то Шз = 1. Значение этого коэффициента будет разъяснено ниже.
2. Асимметричные связи. Переключательная активность
Рассмотрим ансамбль из трёх элементов, объединённых взаимными возбуждающими связями. Матрица коэффициентов связей ансамбля имеет следующий вид:
0 д1 д2 \
а = д2 0 д1
д1 д2 0 /
Рассмотрим сначала случай однонаправленно связанных элементов, то есть примем д\ = 2, д2 = 0. В этом случае матрица описывает асимметричные связи и задаёт порядок активирования элементов 1-2-3. После удаления множителя 2 в каждом уравнении система примет следующий вид:
' р 1 = (Р(р2) + Р(рз))р1(-1 + 2Р(рз) - Р12), < р2 = (Р(рз) + РЫЫ-1 + 2РЫ - р22), (6)
^ р3 = (Р(р1) + Р(р2))рз(-1 + 2Р(р2) - рз2).
Докажем, что если начальные условия принадлежат кубу 0 < р1, р2, рз < 1, то соответствующая фазовая траектория лежит в этой же области. Рассмотрим эволюцию переменной р1, доказательство для двух других аналогично.
Пусть начальные условия таковы, что Р(р2) = Р(рз) = 0. Тогда из первого уравнения системы (6) следует, что = 0 ^ р1 = const. Если Р(р2) = 1, Р(рз) = 0, то р'1 = р1(—1 - р12) и в этом случае р1 монотонно стремится к 0. Если Р(р2) = 0, Р(рз) = 1, то р'1 = р1(1 - р12) и р1 монотонно стремятся к 1. Если Р(р2) = 1, Р(рз) = 1, то р'1 = 2р1(1 - р12) и р1 также монотонно стремятся к 1. Следовательно фазовые траектории действительно лежат в области 0 < р1, р2, рз < 1.
Найдём состояния равновесия системы в области 0 < р1, р2, рз < 1. Для этого рассмотрим все возможные варианты начальных условий.
1. Пусть начальные условия таковы, что Р(Р1) = Р(Р2) = Р(рз) = 0. При этом должны выполняться неравенства 0 < Р1, Р2, Рз < Я. Система примет вид
Р1 = Р2 = Рз = 0. (7)
Отсюда следует, что любая точка (Р1, р2, р3) € {0 < Р1, Р2, Р3 < я} является устойчивым состоянием равновесия.
2. Пусть начальные условия таковы, что Р(Р1) = 1, Р(Р2) = Р(рз) = 0. Эти условия выполнены в области я < р1 < 1, 0 < р2, Рз < Я. Система примет вид
' Р1 = 0,
< р2 = р2(1 - Р22), (8) Рз = Рз( 1 - Рз2)-
Данная система в рассматриваемой области имеет состояние равновесия (С, 0, 0), где я < С < 1. Для точки (С, 0,0) характеристические корни такие: А.1 = 0, А,2 = 1, Аз = -1. Плоскость р1 = С, я < С < 1 - инвариантное множество. В этой плоскости точка (0, 0) будет седлом.
3. Пусть начальные условия таковы, что Р(Р1) = Р(Р2) = 1, Р(рз) = 0. Эти условия выполнены в области я < р1, Р2 < 1, 0 < рз < я. В этом случае система становится такой:
р1 = Р1(-1 - Р12),
< р2 = Р2(1 - Р22), (9)
рз = 2рз(1 - рз2).
Ни одно из её состояния равновесия: (0, 0, 0), (0,1, 0), (0, 0,1), (0,1,1) не принадлежит области я < р1, р2 < 1, 0 < рз < я. Система (9) описывает переходный процесс к случаю, аналогичному предыдущему.
4. Случай Р(Р1) = Р(Р2) = Р(рз) = 1. Эти условия выполняются в области Я < Р1, Р2, Рз < 1. Система примет такой вид
' Р1 = 2р1 (1 - р12),
< р2 = 2р2(1 - Р22), (10) Рз = 2рз(1 - Рз2).
Она имеет одно состояние равновесия (1, 1, 1) с характеристическими числами А,1 = -4, А2 = -4, Аз = -4. То есть это устойчивый узел.
5. Случаи Р(р1) = 0, Р(р2) = 1, Р(рз) =0 и Р(р^ = Р(Р2) = 0, Р(рз) = 1 рассматриваются аналогично случаю 2.
6. Случаи Р(р1) = 1, Р(р2) = 0, Р(рз) = 1 и Р(р1) = 0, Р(Р2) = Р(рз) = 1 рассматриваются аналогично случаю 3.
Из всего сказанного следует, что фазовое пространство системы (6) делится на 3 области.
1) Любая точка из области 0 < р1, р2, рз < Я является неизолированным состоянием равновесия.
2) Любая точка из области д < р1, р2, рз < 1 лежит на траектории, стремящейся к устойчивому узлу (1, 1, 1).
3) Любая другая точка из области 0 < р1, р2, рз < 1 лежит на траектории, движение по которой происходит в плоскостях р1 = С1, р2 = С2, рз = Сз,
д < С1,С2,Сз < 1.
Рассмотрим детально динамику данной системы. Пусть заданы начальные условия (0.5, 0.5, 1)1.
Будем обозначать активное состояние элемента буквой «а», пассивное - буквой «р» (напомним, что г-й элемент считаем активным, если р^ > д). Запись (р, р, а) будет означать, что первый и второй элементы пассивны, третий - активен. Запись (р, р, а) — (а, р, а) будет означать, что в начале процесса был активен третий элемент, а в конце - первый и третий.
1. При указанных начальных условиях система будет иметь вид
' р1 = р1(1 - Pl2), < р2 = р2( 1 - р22), (11) , рз = 0.
Фазовая траектория будет двигаться в инвариантной плоскости рз = 1 к точке (1, 0,1). Отсюда следует, что на этом этапе будет происходить процесс (р, р, а) —
(а, р, а).
2. После того, как первый элемент активируется, система примет вид
' р1 = р1(1 - Pl2), < р2 = 2р2(1 - р22), (12) , рз = рз(-1 - рз2).
Как видно из третьего уравнения, рз устремится к 0. Следовательно, на этом этапе происходит процесс (а, р, а) — (а, р, р)2.
3. Следующий вид системы
' р1 =0,
< р2 = р2(1 - р22), (13)
, рз = рз(-1 - рз2).
На этом этапе фазовая траектория будет двигаться к точке (С, 1, 0) в плоскости р1 = С, д < С < 1. Процесс, происходящий на этом этапе, (а, р, р) — (а, а, р).
1Динамика будет такой же при любых начальных условиях из области 0 < рг < 1, кроме точек 0 < рг < Я и д < р4 < 1.
2Как видно из второго уравнения системы (12), на этом этапе р2 ^ 1. Но время активации много больше времени деактивации, поэтому мы не учитываем этот процесс.
4. Как только второй элемент станет активным, система примет вид
Р1 = Р1(—1 - Р12),
Р2 = Р2(1 — Р22), (14)
Рз = 2рз(1 — Рз2).
Как видно из первого уравнения, р1 устремится к 0. Следовательно, на этом этапе идёт процесс (а, а, р) ^ (р, а, р).
5. На следующем этапе, к которому система придёт после деактивации первого элемента, она примет вид
Р1 = Р1(—1 — Р12),
р2 = 0, (15)
Рз = Рз(1 — Рз2)-
Фазовая траектория будет двигаться в инвариантной плоскости р2 = С к точке (0, С, 1). То есть будет идти процесс (р, а, р) ^ (р, а, а).
6. По окончании активации третьего элемента система придёт к виду
Р1 = 2р1(1 — Р12),
Р2 = Р2 — 1 — Р22), (16)
Рз = Рз(1 — Рз2)-
Из второго уравнения следует, что р2 устремится к 0. Будет идти процесс (р, а, а) ^
^ (р, р, а).
7. Следующим видом системы станет вид из пункта 1.
Таким образом, динамика системы представляет собой процесс циклического изменения активности элементов.
На рис. 1 представлен фазовый портрет для начальных условий из рассмотренного примера, которые обозначены точкой.
На рис. 2 представлены временные реализации для тех же начальных условий. С ростом времени длительности активной фазы элементов увеличиваются. Напомним, что аналогичный процесс происходит в ансамблях асимметрично ингибиторно связанных элементов [13], [1].
На рис. 3 показана зависимость времени активности элементов от порядкового номера активации к для различных значений порога функции Хеви-сайда д.
Рис. 1. Фазовый портрет системы при начальных условиях (0.5, 0.5, 1). Изображающая точка движется поочерёдно в плоскостях р4 = С (д < С < 1), причём С ^ д при t ^ то (если взять начальные условия вида (а,Ь,д), то фазовая траектория будет лежать в плоскостях р4 = д). В этих плоскостях в точках с нулевыми координатами располагаются седловые состояния равновесия с сепаратрисами, направленными по координатным осям. Траектория изображающей точки стремится к этим сепаратрисам
Рг 0.8 0.6 0.4 0.2 0
О
20.0 40.0 60.0 80.0 Рис. 2. Временные реализации колебаний элементов системы. Время активности каждого элемента увеличивается при каждой последующей активации
д/
25.0 20.0 15.0 10.0 5.0 О
_ / / /
/ / / /
/ / ✓ ✓ / / ✓ ✓ ^ - ^ = 0.99 ......... 0.999 — 0.9999 - 0.99999
Такой тип движения обусловливается наличием в фазовом пространстве устойчивого гетероклинического цикла, состоящего из последовательности сед-ловых состояний равновесия и соединяющих их сепаратрис. Каждая траектория, соединяющая сёдла, состоит из двух рёбер куба (см. рис. 1).
Найдём время активности элементов теоретически. Активности г-го элемента соответствует движение по грани куба рг = С (д < С < 1) в фазовом пространстве. Динамика движения по этой грани задаётся системой двух ОДУ Например, движение по грани р1 = д зада-
3
ётся системой
р2 = р2(1 - р22), рз = рз(-1 - рз2)
(17)
0 5.0 10.0 15.0 к Рис. 3. Время активности элемента системы при различных значениях порога функции Хевисайда д. Данные получены путём численного решения
с начальными условиями р2 = а, рз = д при Ь = 0. Движение по грани происходит до точки с координатами р2 = д, рз = Ь при Ь = Т1 (в этот момент происходит переход на другую грань).
Решение системы таково
v/Г—
а2 р2
V1 - р2 а
Уг + рз д
у/1 + д2 рз
Из первого уравнения находим т1
vI-
(18)
от1
а2 д
Из условия рз(т1)
У/г - д2 а
Ь получаем уравнение
VI + Ь2 д
л/ггд2 Ь
В итоге находим зависимость Ь от а
- = е
,Т1
Ь
/
а2(1 - д2)
1 + д2 - 2а2'
(19)
(20)
(21)
3Движение по двум другим граням задаётся такими же системами с точностью до циклической перестановки переменных.
г
е
г
е
Время активности следующего элемента т2 будет задаваться выражением
VI - Ь2 а
еХ2 = , • а- (22)
Тогда выражение для разности между последовательными длительностями активностей будет иметь вид
a
Т2 - Ti = ln :. (23)
Ьл/Г-
a
С учётом (21) эта формула преобразуется к виду
T2-T1 = lV1 + Y^ - гЪ ю - 1 (24)
Данная модель воспроизводит последовательную активность, наблюдаемую в ансамблях, состоящих из биологически релевантных моделей типа Ходжкина-Хаксли [14], и с ее помощью возможно аналитическое изучение этого эффекта.
3. Система с симметричными возбуждающими связями
Рассмотрим теперь систему с такой матрицей коэффициентов связей
0 2 е
G = е 0 2
2 е 0
где 0 < е < 2. Система примет вид
( Рi = (2 + e)(F(Р2) + F(ps))pi(-l + (eF(Р2) + 2F(ps))/Mi - pi2),
J p2 = (2 + e)(F(рз) + F(pi))p2(-1 + (eF(рз) + 2F(pi))/M2 - Р22), (25)
{ p3 = (2 + e)(F(pi) + F(Р2))Рз(-1 + (eF(pi) + 2FЫ)/Мз - Рз2).
Найдём состояния равновесия в области 0 < pi; p2, Рз < 1. Так же, как и ранее, для этого разделим её на несколько частей, в каждой из которых F(p^) = const, i е {1, 2, 3}.
F(pi) = F(p2) = F(рз) = 0. При этом состояния равновесия должны располагаться в кубе 0 < pi; p2, Рз < q. Система примет вид
p'i = Р2 = Рз = 0. (26)
Отсюда следует, что любая точка (pi; p2, Рз) е {0 < pi; p2, Рз < q} является устойчивым состоянием равновесия.
F(pi) = 1, F(p2) = F(рз) = 0. Эти условия выполнены в области q < р1 < 1, 0 < p2, р3 < q. Система преобразуется к виду
pi = 0,
р2 = (2 + е)р2(1 - Р22), (27)
рз = (2 + £)рз(£ - 1 - рз2)-
Эта система в рассматриваемой области имеет состояние равновесия (C, 0, 0), где q < C < 1. При выполнении условия 1 < £ < q2 + 1 = 1.998001 будет ещё одно состояние равновесия (C, 0, л/£ — 1). Для точки (C, 0, 0) характеристические корни такие: = 0, Х2 = 2 + £, Хз = (2 + £)(£ — 1). Для точки (C, 0, л/£ — 1) - Xi = 0, ^2 = 2 + £, Хз = —2(2 + £)(£ — 1) < 0. Таким образом, £ =1 и £ = 1.998001 - бифуркационные значения параметра.
Плоскость р1 = const - инвариантное множество. В этой плоскости точка (0, 0) будет неустойчивым узлом, а точка (0, у/Т—Г) - седлом в случае 1 < £ < q2 + 1. При 0 < £ < 1 точка (0, 0) будет седлом.
F(р1) = F(р2) = 1, F(рз) = 0. Эти условия выполнены в области q < Ръ р2 < 1, 0 < рз < q. В этом случае система становится такой
р1 = (2 + £)р1(£ — 1 — р12),
р2 =(2 + £)Р2(1 — Р22), (28)
Рз = 2(2 + £)рз(£/2 — Рз2)-
При выполнении условия q2 + 1 < £ < 2q2 система будет иметь одно состояние равновесия (V£ — 1, 1, ^/£/2). Его характеристические корни: Х1 = — 2(£ — 1), Х2 = —2, Хз = —£. То есть состояние равновесия является устойчивым узлом. Система (28) описывает переходный процесс к случаю, аналогичному предыдущему.
Случай F(р1) = F(р2) = F(рз) = 1. Эти условия выполняются в области q < р1, р2, рз < 1. Система примет такой вид
р1 = 2(2 + £)р1(£/2 — Р12),
р2 = 2(2 + £)р2(£/2 — Р22), (29)
Рз = 2(2 + £)рз(£/2 — Рз2)-
При выполнении условия у/£/2 > q система будет иметь одно состояние равновесия £/2, у^/2, у^/2). Его характеристические корни: Х1 = —£, Х2 = —£, Хз = —£. Состояние равновесия - устойчивый узел4.
4В этом случае важную роль играет коэффициент Mi: так как все элементы активны, то M1 = M2 = M3 = 2, поэтому устойчивый узел имеет координаты е/2, у/е/2, ^/е/2). При отсутствии коэффициента Mi устойчивый узел находился бы в точке (л/ё, л/ё, л/ё) и при е > 1 изображающая точка вышла бы из кубической области 0 < р1; p2 , Рз < 1.
Случаи Р(р1 )=0, Р(р2)= 1, Р(рз)=0 и Р(р:) = Р(р2) = 0, Р(рз) = 1 рассматриваются аналогично случаю 2.
Случаи Р(р:)= 1, Р(р2)= 0, Р(рз)= 1 и Р(р:) = 0, Р(р2) = Р(рз) = 1 рассматриваются аналогично случаю 3.
Рассмотрим процесс переключательной активности этого ансамбля при условии е — 1 < д. Предположим, что были заданы начальные условия (д, 0.1, 0.1), тогда система примет вид (27). Из второго уравнения этой системы следует, что р2 устремится к единице. Как только р2 достигнет порога функции Хевисайда (напомним, что он у нас равен д = 0.999), система примет вид (28). Как видно из первого уравнения этой системы, начнёт уменьшаться, стремясь к л/е — 1. Но как только станет меньше порогового значения, система примет вид
р 1 = (2 + е)р:(е — 1 — р:2), р 2 = 0, (30)
р з = (2 + е)рз(1 — рз2).
Следовательно, теперь к единице устремится рз, а продолжит движение к значению Vе — 1. Далее, после превышения рз порогового значения, система примет вид
р 1 = 2(2 + е)р1 (е/2 — р:2),
р 2 = (2 + е)р2(е — 1 — р22), (31)
р з = (2 + е)рз(1 — рз2).
Переменная р2 начнёт уменьшаться, стремясь к Vе — 1, и после того, как она станет меньше порогового значения, система примет вид
Г р! = (2 + е)р:(1 — рЛ
! р2 = (2 + е)р2(е — 1 — р22), (32)
I рз = 0.
Как видно из первого уравнения, устремится к единице. В следующий раз система сменит вид после того, как станет равно пороговому значению
( р! = (2 + е)р:(е — 1 — р^),
\ р2 = 2(2 + е)р2(е/2 — р22), (33)
I рз = (2 + е)рз(е — 1 — рз2).
Переменная рз станет уменьшаться. После того, как она станет меньше порогового значения, система придёт к исходному виду, и процесс повторится.
Из вышесказанного следует, что колебания элементов данного ансамбля будут происходить в диапазоне (Vе — 1, д]. На рис. 4 приведены временные реализации колебаний элементов системы при различных значениях е.
Рассмотрим теперь ансамбль при условии д2 + 1 < е < 2. Если активных элементов нет, то есть Р(Р1) = Р(р2) = Р(рз) = 0, то система примет вид р 1 = р2 = р3 = 0.
Если активен только один элемент, например, первый (то есть Р(р1) = 1, Р(р2) = Р(р3) = 0), то система будет выглядеть так
р1 = 0,
р2 = (2 + е)р2(1 - р22), (34)
рз = (2 + е)рз(е - 1 - рз2).
В этом случае через некоторое время станут активны все 3 элемента.
Если активны любые 2 элемента, например, первый и третий (то есть Р(р1) = 1, Р(р2) = 0, Р(р3) = 1), то система преобразуется к виду
р1 = (2 + е)р1(1 - р12),
р2 = 2(2 + е)р2(е/2 - р22), (35)
рз = (2 + е)рз(е - 1 - рз2).
Так как у/е/2 > е - 1 > д, то в этом случае также станут активны все 3 элемента.
Если активны все 3 элемента (Р(р1) = Р(р2) = Р(рз) = 1), система примет вид
р1 = 2(2 + е)р1(е/2 - р12),
р2 = 2(2 + е)р2(е/2 - р22), (36)
рз = 2(2 + е)рз(е/2 - рз2).
В этом случае изображающая точка будет стремиться к устойчивому узлу е/2, ^е/2, ^е/2) и все элементы останутся активными.
Из вышесказанного видно, что при таких значениях параметра е у системы есть 2 устойчивых состояния: когда все элементы неактивны и когда все активны (р1 = р2 = рз = л/е/2). Последовательной передачи активности здесь нет.
Найдём возможные типы движения в фазовом пространстве в завимо-сти от значения е (0 < е < д2 + 1). Зададим начальные условия в точке (д, а, д)5 и построим точечное отображение ребра куба р1 = рз = д, 0 < р2 < д в ребро р1 = р2 = д, 0 < рз < д. Динамика системы в инвариантной плоскости р1 = д задаётся системой
( р2 = р2(1 - р22),
2 (37)
I рз = рз(е - 1 - рз2).
5 Будем рассматривать начальные условия только такого вида, так как при любых начальных условиях, при которых активен только один элемент, координаты изображающей точки через некоторое время примут такой вид.
Рис. 4. Временные реализации колебаний системы при различных значениях е: а -1.5, б -1.9, в - 1.98. При увеличении е амплитуда колебаний уменьшается до нуля
При е = 1 её решение, удовлетворяющее начальным условиям, таково
P2V/Í—
2
«V1 - Р22 1 1
з"
(38)
2р
г.
Конечная точка - это точка с координатами р2 = д, рз = Ь6. Зависимость Ь(а) задаётся формулой
Ь
\
1
1 д2(1 - а2)-
(39)
д2
а2(1 - д2)
График этой функции имеет 3 точки пересечения с прямой Ь = а: на концах отрезка [0, д] и во внутренней точке этого отрезка, которой соответствует устойчивый предельный цикл7.
При е =1 решение системы будет таким
' Р2л/1 - а2 .
—/ = е ■ аV 1 - Р22
РзУ |д2 - е +1 , дл/|рз2 - е +1|
(40)
,(£-1)1
6а, Ь - координаты точек пересечения фазовой траектории с рёбрами куба.
7Для вывода о наличии или отсутствии предельного цикла необязательно рассматривать все 3 инвариантные плоскости = д, так как в двух других динамика будет в точности такая же.
г
е
для Ь
Отсюда получаем выражение
Ь
\
д2(е -1)
д2-\д2 - е+1|
(а2(1 - д2) \ \д2(1 - а2))
£- 1 •
Рис. 5. График зависимости минимального уровня спайковой активности от е. Пунктиром показан график функции Vе — 1
д2(1 - а2)
' (41)
Функция Ь(а) возрастающая. Её график пересекается с осью ординат в точке Vе - 1 при 1 < е < д2 + 1 и в 0 при 0 < е < 1. Он также проходит через точку (д, д) и имеет ещё одну точку пересечения (а*,Ь*) с прямой Ь = а (то есть а* = Ь*) при 0 < а* < д. Этой точке соответствует устойчивый предельный цикл.
Таким образом, в фазовом пространстве системы при 0 < е < < д2 + 1 существует устойчивый предельный цикл, в окрестности которого будет двигаться фазовая точка. В этом случае Ь* будет минимальным уровнем спайковой активности. На рис. 5 приведен график зависимости Ь* от е.
Время активности элементов можно найти из первой формулы системы (40) и формулы (41). Для этого нужно задать начальное значение а0 и воспользоваться этими формулами
Рис. 6. Графики зависимости времени активности элементов от порядкового номера активации для различных значений е
М
г+1
аг+1
1п
gv/Г-07
\
д (е -1)
(42)
д2 - \д2 - е + 1\
(а2 (1 - д2) N \д2(1 - а^2) /
£- 1'
На рис. 6 представлены графики зависимости времени активности элементов от порядкового номера активации для различных значений е.
2
а
д
4. Влияние на систему малого параметра
В реальных ситуациях на поведение нейронных ансамблей влияют внешние возмущения. Рассмотрим влияние малого параметра ^ на динамику исходной системы с асимметричными связями.
Рис. 8. Временные реализации колебаний системы при наличии малого возмущающего паРис. 7. Фазовый портрет системы при нали- раметра. Время активности элементов не уве-чии малого возмущающего параметра. В этом личивается. Такое движение соответствует дви-случае траектория изображающей точки пере- жению в окрестности устойчивого предельного стаёт стремиться к рёбрам куба цикла
Возмущённая система имеет вид
Р1 = (Р(р2) + Р(рз))(р1(-1 + 2Р(рз) - Р12) + и),
Р2 = (Р(Рз) + Р(Р1 ))(р2(-1 + 2Р(Р1) - Р22) + и), (43)
Р3 = (Р(Р1) + Р(Р2))(рз(-1 + 2Р(Р2) - Рз2) + и)-
Из проведённого численного расчёта видно, что в этом случае нет увеличения времени, которое фазовая точка проводит в окрестности сёдел. Фазовый портрет при и = 0-001 приведён на рис. 7. На рис. 8 приведены соответствующие временные реализации.
Таким образом, при введении в систему малого параметра и происходит разрушение гетероклинического контура и появляется устойчивый предельный цикл.
Выводы
В данной статье была предложена новая модель нейронного ансамбля. Математическая модель элементов этого ансамбля задаётся обобщённым уравнением Лотки-Вольтерры и элементы возбуждающе связаны друг с другом.
Рассмотрены модели с симметричными и несимметричными связями. Установлено, что при несимметричных связях в фазовом пространстве системы существует устойчивый гетероклинический цикл. Устойчивый гете-роклинический цикл также существует в ансамбле ингибиторно связанных нейронов. При симметричных связях происходит разрушение гетероклинического цикла и появляется устойчивый предельный цикл.
Кроме того, проведено исследование возмущённой системы. Выяснилось, что при введении малого возмущающего параметра происходит разрушение устойчивого гетероклинического цикла и появляется устойчивый предельный цикл.
Несмотря на то, что исследованная система является феноменологической моделью нейронного ансамбля, она демонстрирует основные эффекты, которые наблюдаются в ансамблях, состоящих из биологически релевантных моделей типа Ходжкина-Хаксли.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» (соглашение № 14.В37.21.0863).
Библиографический список
1. Rabinovich M.I., VaronaP., Selverston A.I., and Abarbanel H.D.I. Dynamical principles in neuroscience // Review of modern physics. 2006. Vol. 78. 1213.
2. 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.
3. 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 // Neural comput. 2004. Vol. 16. P. 999.
4. LeviR., VaronaP., Arshavsky Y.I., Rabinovich M.I., and Selverstone A.I. Dual sensory-motor function for a molluskan statocyst network // Neurophysiol J. 2004. Vol. 91. P. 336.
5. Guckenheimer John, and Holmes Philip. Structurally stable heteroclinic cycles // Math. Proc. Camb. Phil. Soc. 1988. Vol. 103. P. 18.
6. Emily Stone, and Philip Holmes. Random perturbations of heteroclinic attrac-tors // SIAM Journal on Applied Mathematics. 1990. Vol. 50, № 3. P. 726.
7. Postlethwaite Claire M., and Dawes Jonathan H.P. Resonance bifurcations from robust homoclinic cycles // Nonlinearity. 2010. Vol. 23. P. 621.
8. Driesse Ramon, and Ale Jan Homburg. Resonance bifurcation from homo-clinic cycles // Differential Equations. 2009. Vol. 246. P. 2681.
9. Seliger P., Tsimring L.S., and Rabinovich M.I. Dynamics-based sequential memory: Winnerless competition of patterns // Physical Review E. 2003. Vol. 67. 011905.
10. Afraimovich VS., Rabinovich M.I., and Varona P. Heteroclinic contours in neural ensembles and the winnerless competition principle // Int. J. of Bifurcation and Chaos. 2004. Vol. 14, № 4.
11. Afraimovich VS., Zhigulin V.P., and Rabinovich M.I. On the origin of reproducible sequential activity in neural circuits // Chaos. 2004. Vol. 14, № 4.
12. Rabinovich M., Volkovskii A., Lecanda P., Huerta R., Abarbanel H.D.I., and Laurent G. Dynamical encoding by networks of competing neuron groups: Winnerless competition // Physical review letters. 2001. Vol. 87, № 6.
13. Рабинович М.И., Мюезинолу М.К. Нелинейная динамика мозга: Эмоции и интеллектуальная деятельность // Успехи физических наук. 2010. Vol. 180, № 4.
14. Komarov M.A., Osipov G.V. and Suykens J.A.K. Sequentially activated groups in neural networks // EPL. 2009. Vol. 86. 60006.
Нижегородский госуниверситет Поступила в редакцию 23.04.2013
им. Н.И. Лобачевского После доработки 9.07.2013
SEQUENTIAL ACTIVITY IN NEURONAL ENSEMBLES WITH EXCITATORY COUPLINGS
A. G. Korotkov, G. V. Osipov
A new model of neurons like elements is suggested in the paper. The model is based on the generalized Lottka-Volterra model with excitatory coupling. The study is motivated by the fact that the excitatory couplings are the dominating type of interactions between neurons in the human brain. It is shown in the paper that there are two regimes exist in such ensemble of oscillators in dependence on the coupling between the elements: the regime with stable heteroclinical cycle and the regime with stable limit cycle.
Keywords: Neuronal ensembles with excitatory coupling, Lottka-Volterra model, numerical simulation.
Коротков Александр Геннадьевич - родился в Нижнем Новгороде (1981). Студент 6 курса вечернего отделения факультета вычислительной математики и кибернетики Нижегородского государственного университета.
603950 Россия, Нижний Новгород, пр. Гагарина, 23 Нижегородский государственный университет им. Н.И. Лобачевского E-mail: [email protected]
