Нелинейные динамические модели нейронов: обзор Текст научной статьи по специальности «Математика»

УДК 621.373.1

Нелинейные динамические модели нейронов: Обзор

А. С. Дмитричев, Д. В. Касаткин, В. В. Клиньшов, С. Ю. Кириллов, О. В. Масленников, Д. С. Щапин, В. И. Некоркин

Институт прикладной физики РАН Россия, 603950 Нижний Новгород, БОКС-120, Ульянова, 46 E-mail: admitry@neuron.appl.sci-nnov.ru, kasatkin@appl.sci-nnov.ru, vladimir.klinshov@gmail.com, skirillov@ipfran.ru, olmaov@ipfran.ru, shapinds@mail.ru, vnekorkin@neuron.appl.sci-nnov.ru Поступила в редакцию 01.06.2018, принята к публикации 28.06.2018

Тема исследования. Представлен обзор основных динамических моделей нейронной активности и обсуждаются индивидуальные особенности их поведения, которые могут быть в последующем использованы как основа при разработке и построении различных конфигураций нейронных сетей. Работа содержит как новые оригинальные результаты, так и обобщение уже известных, опубликованных ранее в разных журналах. Цель -познакомить читателя с базовыми динамическими свойствами нейронов, такими как наличие состояния покоя и генерация потенциала действия; сформировать у него общее представление о динамических механизмах, лежащих в основе отмеченных свойств и используемых при построении моделей нейронной активности различного уровня детализации. Исследуемые модели. С математической точки зрения модели нейронов делятся на два класса. Первый класс представлен моделями с непрерывным временем в форме систем обыкновенных дифференциальных уравнений. Моделям этого класса посвящена вторая часть данного обзора. Открывает эту часть наиболее детализированная модель Ходжкина-Хаксли, являющаяся канонической моделью нейронной активности в нелинейной динамике. Затем приводятся упрощенные модели - двумерная модель Моррис-Лекара для спайкинга и трехмерная модель Хиндмарша-Роуза для бёрстинга. Наиболее подробно описана модель ФитцХью-Нагумо, для которой проведен детальный бифуркационный анализ. Также приведены модели, описывающие нейроны со специфическими свойствами - модель нейрона с постдеполяризацией и модель нейрона нижних олив. Завершает эту часть наиболее простая модель типа «накопление-сброс». Второй класс образуют модели с дискретным временем, представляющие собой точечные отображения. Такие модели в последнее время приобретают все большую популярность в виду богатства демонстрируемых динамических режимов и простоты численного моделирования. Моделям этого класса посвящена третья часть данного обзора. В частности, приведены такие модели как модель Киалво, модель Ижикевича, модель Рулькова и модель Курбажа-Некоркина. Результаты. Изложены базовые физические принципы построения математических моделей нейронной активности, основанные на ионном транспорте. На примере модели ФитцХью-Нагумо изучены основные свойства и механизмы возникновения режимов мультипорогового возбуждения в нейронах. Раскрыт механизм формирования бёрстовых колебаний в модели Хиндмарша-Роуза. Описан динамический механизм временного понижения порога возбуждения и возникновения периодических колебаний в модели нейрона с постдеполяризацией. Описано формирование в нейронах нижних олив (Ca2+)- и (№2+)-зависимых спайков. Описаны динамические механизмы

формирования основных регулярных и хаотических режимов нейронной активности в дискретных моделях Киалво, Ижикевича, Рулькова и Курбажа-Некоркина. Обсуждение. В Заключении кратко резюмируется содержание обзора.

Ключевые слова: динамические системы, нейроны, бифуркации.

https://doi.org/10.18500/0869-6632-2018-26-4-5-58

Образец цитирования: Дмитричев А.С., Касаткин Д.В., Клиньшов В.В., Кириллов С.Ю., Масленников О.В., Щапин Д.С., Некоркин В.И. Нелинейные динамические модели нейронов // Известия вузов. Прикладная нелинейная динамика. 2018. Т. 26, № 4. С. 5-58. https://doi.org/10.18500/0869-6632-2018-226-4-5-58

Nonlinear dynamical models of neurons: Review

A. S. Dmitrichev, D. V. ^sat^n, V. V. Klinshov, S. Yu. Kirillov, O. V. Maslennikov, D. S. Shchapin, V. I. Nekorkin

Institute of Applied Physics of RAS 46, Ul'yanov st., 603950 Nizhny Novgorod, Russia E-mail: admitry@neuron.appl.sci-nnov.ru, kasatkin@appl.sci-nnov.ru, vladimir.klinshov@gmail.com, skirillov@ipfran.ru, olmaov@ipfran.ru, shapinds@mail.ru, vnekorkin@neuron.appl.sci-nnov.ru Received 01.06.2018, accepted for publication 28.06.2018

Topic. A review of the basic dynamical models of neural activity is presented and individual features of their behavior are discussed, which can be used as a basis for the subsequent development and construction of various configurations of neural networks. The work contains both new original results and generalization of already known ones published earlier in different journals. Aim is to familiarize the reader with the basic dynamical properties of neurons, such as the existence of a rest state and the generation of the action potential; to outline the dynamical mechanisms underlying these properties which are used in the development of neural models with various levels of detailing. Investigated models. From the mathematical point of view, neuron models are divided into two classes. The first class is represented by models with continuous time described by ordinary differential equations. The section devoted to continuous-time models starts from the most detailed Hodgkin-Huxley model, which is a canonical model for neural activity in nonlinear dynamics. Further we describe simplified models, such as a two-dimensional model of Morris-Lecar for spiking and a three-dimensional model of Hindmarsh-Rose for bursting. The FitzHugh-Nagumo model is described in detail, and detailed bifurcation analysis is presented. We also present models for neurons with specific properties, namely a neuron with afterdepolarization and an inferior olives neuron. The last and the simplest model is the «integrate-fire» model. The second class of neural models are systems with discrete time represented by discrete maps. Such models have recently gained increasing popularity due to the richness of the demonstrated dynamics and the ease of numerical simulations. We describe such models as the Chialvo model, the Izhikevich model, the Rulkov model, and the Courbage-Nekorkin model. Results. The basic physical principles underlying the construction of mathematical models of neural activity, based on ion transport, are outlined. Using the FitzHugh-Nagumo model as an example, the main properties and mechanisms of the emergence of multithreshold excitation regimes in neurons are described. The mechanism of formation of burst oscillations in the Hindmarsh-Rose model is outlined. A dynamic mechanism for temporal decline of the excitation threshold and the emergence of periodic oscillations in a neuron with afterdepolarization are described. The formation of (Ca2+)- and (Na2+)-dependent spikes in inferior olive neurons is described. Dynamic mechanisms of formation of the major regular and chaotic regimes of neural activity in discrete models of Chialvo, Izhikevich, Rulkov and Courbage-Nekorkin are described. Discussion. In the Conclusion we briefly summarize the content of the survey.

Keywords: dynamical systems, neurons, bifurcations.

https://doi.org/10.18500/0869-6632-2018-26-4-5-58

References: Dmitrichev A.S., Kasatkin D.V., Klinshov V.V., Kirillov S.Yu., Maslennikov O.V., Shchapin D.S., Nekorkin V.I. Nonlinear dynamical models of neurons. Izvestiya VUZ, Applied Nonlinear Dynamics, 2018, vol. 26, iss. 4, pp. 5-58. https://doi.org/ 10.18500/0869-6632-201826-4-5-58

Введение

Нейроны (или нервные клетки) являются основными структурными и функциональными единицами нервной системы. Как и любая биологическая клетка, нейрон представляет собой клеточное тело размером порядка 10-100 рш, ограниченное снаружи мембраной - непроницаемым двойным слоем липидов. Отличительным свойством мембраны нервных клеток является наличие разности электрических потенциалов между внутренней и наружной поверхностями. Эта разность, называемая мембранным потенциалом, имеет в состояния покоя некоторое постоянное отрицательное значение и может изменяться во времени. Генерируемые мембраной электрические импульсы служат для кодирования и передачи информации в нейронных сетях. Передача между нейронами электрических сигналов связана с их распространением по специальным отросткам, отходящим то тела клетки в виде тонких трубчатых волокон. Известно два вида таких отростков - аксон и дендриты. Аксон -длинный отросток, предназначенный для передачи электрических сигналов от тела нейрона к другим клеткам. Дендриты - обычно более короткие, разветвленные отростки, через которые на нейрон поступают входящие сигналы. Аксон заканчивается многочисленными терминальными волокнами, которые могут контактировать с дендритами, аксоном или телом другого нейрона. Этот контакт осуществляют специальные устройства - синапсы (рис. 1). Синапс состоит из двух мембран (преси-наптической и постсинаптической) контактирующих нейронов. Основными видами синапсов являются так называемые химические и электрические синапсы. В электрическом синапсе расстояние между мембранами пресинаптического и постсинап-тического нейронов составляет порядка 2-4 пш, а их взаимодействие осуществляется с помощью особых белковых образований. Совершенно иначе происходит контакт между нейронами в случае химического синапса. Мембраны взаимодействующих нейронов в этом случае удалены друг от друга на расстояние порядка 20-40 пш и образуют так называемую синаптическую щель. Когда электрический сигнал достигает пресинаптического окончания из его мембраны высвобождается особое вещество - нейромедиатор. Нейромедиатор поступает в синаптическую щель, пройдя которую начинает действовать на рецепторы постсинаптической мембраны, изменяя её потенциал, и тем самым осуществляет связь между нейронами. Синапсы различаются как по типу используемого нейромедиатора, так и по функциональным свойствам. Если под действием нейромедиатора мембранный потенциал постсинап-тического нейрона приближается к порогу, то есть происходит деполяризация его мембраны, то синапс называется возбуждающим. В противном случае, когда мембранный потенциал постсинаптического нейрона удаляется от порога и происходит гиперполяризация его мембраны, синапс называется тормозным или подавляющим.

Рис. 1. Схематическое изображение нейрона и синапса Fig. 1. Schematic representation of a neuron and a synapse

Основной формой электрической активности нейрона является так называемый потенциал действия - электрический импульс или спайк. В невозбужденном состоянии мембрана нейрона поляризована - её внутренняя поверхность имеет отрицательный потенциал по отношению к наружной. Это так называемый потенциал покоя, равный обычно (—70)-(—90) mV. При возникновении потенциала действия сначала наблюдается фаза нарастания, во время которой мембранный потенциал достаточно резко нарастает до пикового значения (+20)-(+40) mV, после чего некоторое время наблюдается пик возбуждения. Затем наступает фаза реполяризации, когда мембраный потенциал вновь снижается до значения, чуть меньшего, чем исходное. Суммарная длительность потенциала действия в нервных клетках составляет порядка 1 ms, после чего ещё в течение 1-3 ms в мембране наблюдаются некоторые остаточные явления. Это так называемый период рефрактерности, во время которого мембрана невозбудима. Характеристики потенциалов действия практически идентичны во всех нервных клетках организма и являются универсальными носителями информации в нейронной сети.

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

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

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

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

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

Впервые понятие формального нейрона было сформулировано в статье Мак-Каллока и Питса [1]. Формальные нейроны представляют собой пороговые элементы, осуществляющие преобразование вход-выход. Такой нейрон получает на входе вектор входных сигналов, предварительно умноженный на вектор весовых коэффициентов, имитирующих действие синапсов. Затем происходит суммирование координат входного вектора и над полученной суммой выполняется нелинейное преобразование. В результате на выходе формального нейрона получается сигнал, который может быть как скалярной переменной, так и вектором. Далее этот сигнал посылается на входы других нейронов, формируя тем самым искусственную нейронную сеть. В настоящее время такие сети широко используются для решения научных и инженерных задач. Однако для моделирования процессов активности мозга, возникающих, в частности, при выполнении когнитивных функций, требуются модели нейронов более адекватные с нейрофизиологической точки зрения. Мы придём к таким моделям, например, если учтём наличие у нейронов ионных каналов, благодаря которым возникают ионные токи и происходит изменение мембранного потенциала. Поскольку эти процессы являются эволюционными, соответствующие им математические модели имеют форму нелинейных динамических систем. В данном обзоре рассматриваются нейронные модели именно этого класса, начиная с известной модели Ходжкина-Хаксли.

Статья организована следующим образом. В Разделе 1 описаны физические принципы, лежащие в основе возбудимости нервных клеток. Раздел 2 посвящен моделям с непрерывным временем. В нем приведены основные сведения о модели Ходжкина-Хаксли. Далее рассмотрены модели Моррис-Лекара, ФитцХью-Нагумо, Хиндмарша-Роуза, модели нейрона с последеполяризацией, модели нейрона

нижних олив и модели типа «накопление-сброс». В Разделе 3 изложены подходы к описанию динамики нейроннов с использованием систем с дискретным временем. Рассмотрены различные моды нейронной активности, воспроизводимые в модели Киалво, модели Ижикевича, модели Рулькова и модели Курбажа-Некоркина. В Заключении кратко обсуждаются математические аспекты представленных в обзоре моделей.

1. Физические основы нейронной возбудимости

Физико-химическую основу возбуждения нервных клеток составляют процессы ионного транспорта через клеточную мембрану. В данном транспорте участвуют в основном четыре вида ионов: КСа2+, С1-. Несмотря на то, что мембрана клетки представляет собой непроницаемый для этих ионов липидный бислой, в нее встроены специализированные белковые системы, которые обеспечивают транспорт ионов через клеточную мембрану. Существует два типа таких систем - ионные насосы и ионные каналы. Ионные насосы (помпы) - интегральные белки, которые обеспечивают активный перенос ионов против градиента концентрации. В результате работы ионных насосов создаются и поддерживаются трансмембранные ионные градиенты: концентрация натрия, кальция и хлора внутри клетки ниже, чем снаружи (в межклеточной жидкости); концентрация калия внутри клетки выше, чем снаружи. Разность ионных концентраций внутри и снаружи клетки приводит к возникновению ненулевых потенциалов реверсии для данных типов нейронов и, следовательно, отличного от нуля мембранного потенциала.

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

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

Свойство клеточной мембраны накапливать на своей внутренней поверхности заряд позволяет рассматривать её как электрический конденсатор определённой емкости. Ионные токи, протекающие через мембрану, приводят к зарядке или разрядке данного конденсатора. Исходя из сказанного, нейрон может быть представлен эквивалентной электрической схемой, изображенной на рис. 2. Записывая уравнение Кирхгоффа

для проходящих через поверхность клеточной мембраны токов, получим следующее уравнение:

Рис. 2. Эквивалентная схема нейронной мембраны

Fig. 2. Electronic equivalent circuit for a neural membrane

N

C— = Iext 1

(1)

j=1

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

(2)

Ij = 9i(V - Ej).

Уравнения (1)-(2) являются основными уравнениями трансмембранного ионного транспорта. Существует широкий класс моделей нейронов, основанных на данных уравнениях. Эти модели отличаются количеством и типом учитываемых ионных каналов, а также способами описания проводимости селективных каналов, которая динамическим образом изменяется во времени благодаря их открытию и закрытию.

2. Модели нейронной активности с непрерывным временем

2.1. Модель Ходжкина-Хаксли. Исторически первой и наиболее известной моделью, основанной на ионном транспорте, является модель, предложенная в 1952 году английскими нейрофизиологами А. Ходжкиным и Э. Хаксли [2]. В своих экспериментах, выполненных на гигантском аксоне кальмара, они исследовали механизмы возникновения и передачи нервного возбуждения. В своих работах Ходжкин и Хаксли показали, что проводимость селективных каналов определяется состоянием их так называемых активирующих и инактивирующих ворот. Активирующие ворота постепенно открываются при повышении мембранного потенциала, а инак-тивирующие ворота, наоборот, при этом закрываются. Проводимость каналов ]-го типа определяется следующим выражением:

д3 = дз ар/ Щ, (3)

где дщ - максимальная проводимость каналов в открытом состоянии; ащ и Ьщ - соответственно активационная и инактивационная переменные, характеризующие вероятность быть в открытом состоянии активирующих и инактивирующих ворот; рщ и

qj - положительные целые числа, определяющие число активирующих и инактиви-рующих ворот на ионный канал. По предположению Ходжкина и Хаксли, динамика активационной и инактивационной переменных описываются дифференциальными уравнениями первого порядка

dx dt

x^(V) - x Тж (V) •

(4)

Здесь х = aj,bj - вектор переменных; х^(У) и х^,(У) определяют стационарные состояния переменных и время релаксации процессов, которые в общем случае зависят от величины мембранного потенциала. Ходжкин и Хаксли установили, что главную роль в процессе ионного транспорта играют ионы натрия и калия. Именно натриевые и калиевые каналы детально описываются в предложенной ими модели, а вклад остальных ионов характеризуется общим током утечки. Ходжкин и Хаксли экспериментально показали, что для натриевых каналов характерно наличие трех активирующих и одних инактивирующих ворот (рма = 3, qNa = 1), а для калиевых -наличие четырех активирующих ворот на канал (рк = 4, qк = 0). Данные числа подбираются эмпирически так, чтобы наилучшим образом описать кинетику токов. Функции х^(У) и х^(У) также подбираются эмпирически. В результате модель Ходжкина-Хаксли принимает следующий вид:

1У

с— = Тех* - дкп4(У - Ук) - дмат3Ь(У - Ума) - дь(У - Уь), (5)

dt dx

ddt = ax(V)(1 - x) - |3x(V)x.

(6)

Здесь У обозначает мембранный потенциал в шУ, отсчитываемый от потенциала покоя; вектор х = п,т, Н, где п и т - активационные переменные для калиевых и натриевых каналов, Н - инактивационная переменная натриевых каналов.

Численное интегрирование системы (5)-(6) при биологически релевантных параметрах показывает, что данная модель адекватно воспроизводит динамику мембранного потенциала нейрона. Основным свойством системы является генерация потенциала действия, происходящая по закону «все или ничего». На рис. 3 приведены

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

Рис. 3. Сравнение экспериментально наблюдаемого потенциала действия в гигантском аксоне кальмара и результатов численного интегрирования модели Ходжкина-Хаксли. Модифицировано из [2]

Fig. 3. Comparison of the action potential experimentally observed in the giant squid axon and the results of numerical integration of the Hodgkin-Huxley model. Modified from [2]

Принципы, положенные в основу модели Ходжкина-Хаксли, были использованы другими авторами для создания целого ряда моделей подобного типа (см., например, [3-5]). Эти модели отличаются количеством и типом учитываемых ионных каналов. Несомненным достоинством моделей типа Ходжкина-Хаксли является детальное описание динамики ионных каналов, биологическая адекватность и наличие соответствия модельных и физиологических параметров. Однако сложность аналитического и численного исследования привела к необходимости использования других, более простых моделей нейронной активности.

2.2. Модель Моррис-Лекара. Другая широко известная модель, опирающаяся на формализм Ходжкина-Хаксли, была предложена К. Моррис и Г. Лека-ром [6] для описания различных режимов осцилляторной активности, наблюдаемых в экспериментах с гигантским мышечным волокном усоногого рака. На основе экспериментальных данных [7] была построена трехмерная динамическая модель, учитывающая наличие двух потенциал-зависимых ионных каналов - кальциевого и калиевого, которые никогда полностью не инактивируются. В работе [6] представлен анализ возможных режимов, предсказываемых теоретической моделью, и их сопоставление с данными, полученными экспериментально. При проведении анализа режимов, отвечающих осцилляторному поведению потенциала мышечного волокна было сделано предположение, что проводимости кальциевого и калиевого каналов имеют различные времена релаксации. Это позволило получить редуцированную модель второго порядка при условии, что проводимость кальциевого канала мгновенно принимает стационарное значение,

(IV

С— = !ехь - дь(V - Vь) - доаМ^ФXV - Vca) - дкп(V - VК), (7) м

(п

Пг

.(V) - п хп(У) -

Переменные V и п имеют тот же смысл, что и в модели Ходжкина-Хаксли. Слагаемые в правой части первого уравнения (7) определяют, соответственно, внешний ток, ток утечки, кальциевый и калиевый токи, текущие через мембрану. Функции М^), Пж^) и тn(V) имеют следующий вид:

М^) = 0.5

п^) = 0.5

1 + tanh

1 + tanh

V - VI V2

V-Уз V4

Тп^) =

ф ео8Ь(

-1

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

Динамика модели (7) была рассмотрена в ряде работ [8-12]. Основное внимание в данных работах преимущественно было сосредоточено на исследовании влияния внешнего тока на характер рождения колебаний. Наиболее подробно исследование бифуркационной структуры модели М-Л в плоскостях различных параметров представлено в [11]. Проведённый в этой работе анализ двухпараметрических бифуркационных диаграмм, где в качестве контрольных параметров выбраны внешний ток Тех1 и одна из величин дса, ф, Уз или У4, позволил выделить области с различным поведением мембранного потенциала и определить сценарии, приводящие к установлению этих режимов. Система (7) может демонстрировать режим периодической генерации импульсов, а также возбудимые режимы различной сложности. В общем случае, поведение системы определяется числом аттракторов, взаимным расположением сепаратрис седла и неустойчивых предельных циклов, существующих на фазовой плоскости при выбранных значениях параметров. В модели М-Л существует два вида аттракторов - устойчивые состояния равновесия и предельные циклы, которые при описании нейронной активности отвечают состоянию покоя и периодической генерации импульсов нейрона. В зависимости от параметров динамика модели М-Л может быть моностабильной, бистабильной и даже тристабильной. Установлено, что система (7) может иметь от одного до трех состояний равновесия, что определяется числом точек пересечения нульклин системы (рис. 4). В случае, когда система имеет три состояния равновесия, одно из них (О2) является седлом, устойчивые сепаратрисы которого задают порог возбуждения нейрона. Простейший пример возбудимого режима представлен на рис. 4, a. В начальный момент времени система находится в состоянии покоя, которое на фазовой плоскости определяется устойчивым состоянием равновесия 0\ и которое является единственным аттрактором системы. При воздействии на систему внешнего стимула происходит смещение изображающей точки на фазовой плоскости из состояния равновесия. После окончания действия стимула система возвращается в состояние покоя. При этом если внешнее воздействие было достаточно сильным, то отклик характеризуется значительным нарастающим изменением переменной У, то есть система генерирует импульс возбуждения. Система (7) может демонстрировать и более сложное возбудимое поведение, когда в ответ на внешнее воздействие наблюдается генерация серии импульсов возбуждения (рис. 4, б). Такое поведение наблюдается при приближении параметров системы к бифуркационным значениям, отвечающим рождению двукратного предельного цикла. Возникновение сложного отклика обусловлено образованием «слоений» на фазовой плоскости, связанным со сгущением фазовых траекторий, являющимся следствием динамической памяти системы о недавнем существовании двукратного предельного цикла. При попадании в одно из таких «слоений», возвращаясь в состояние равновесия О1, траектория системы совершает некоторое число оборотов, что означает генерацию серии импульсов возбуждения.

Режим периодической генерации импульсов возбуждения определяется появлением на фазовой плоскости системы устойчивого предельного цикла (рис. 4 c, й). Режим периодических колебаний может возникать в результате различных бифуркаций [11]: Андронова-Хопфа, двукратного предельного цикла, бифуркации петли сепаратрис седло-узла.

■60 -40 -20 0 20 »,mV 0 20 40 60 i,ms

■60 -40 -20 0 20 f,mV 0 10 20 jgms

Рис. 4. Фазовые портреты и временные реализации модели М-Л при C = 20 ^F/cm2, gL = 0.2 mS/cm2, gca = 4.4 mS/cm2, дк = 8.0 mS/cm2, Vi = 50 mV, Vca = 100 mV, Vk = 70 mV, Vi = -1 mV, V2 = 15 mV, V3 = 2 mV, V4 = 30 mV, ф = 0.05 s_1 и при различных значениях внешнего тока Iext |j,A/cm2: a - 20, b - 27.54842, с -30, d - 40

Fig. 4. Phase portraits and waveform of membrane potential in Morris-Lecar model for C = 20 ^F/cm2, gL = 0.2 mS/cm2, gca = 4.4 mS/cm2, дк = 8.0 mS/cm2, Vi = 50 mV, Vca = 100 mV, Vk = 70 mV, Vi = -1 mV, V2 = 15 mV, V3 = 2 mV, V4 = 30 mV, ф = 0.05 s^1 and for different value of external current Iext и-A/cm2: a - 20, b - 27.54842, с -30, d - 40

Система (7) может демонстрировать бистабильное поведение, одновременно обладая возбудимыми и колебательными свойствами. Пример фазового портрета, отвечающего такому поведению, представлен на рис. 4, e. На фазовой плоскости существуют пара предельных циклов - устойчивый и неустойчивый и три состояния равновесия, одно из которых (О1) также является устойчивым. Если в результате действия стимула изображающая точка оказывается ниже устойчивых сепаратрис седла О2, оставаясь при этом внутри области, ограниченной неустойчивым предельным циклом, то система генерирует импульс возбуждения (см. рис. 4, e). Если параметры стимула таковы, что изображающая точка на фазовой плоскости после его действия оказывается вне области, выделяемой неустойчивым предельным циклом, то система переходит в режим периодических колебаний, отвечающих устойчивому предельному циклу.

Наконец, существуют области параметров, где система имеет три аттрактора [11]: два состояния равновесия и предельный цикл. В этом случае благодаря внешней стимуляции можно вызвать, например, генерацию серий импульсов возбуждения, заканчивающуюся возвращением обратно в состояние покоя, или перевести в режим периодических колебаний.

2.3. Модель ФитцХью-Нагумо. Как было сказано выше, динамика селективного ионного канала характеризуется следующими ключевыми чертами. Во-первых, при фиксированном трансмембранном потенциале V значение тока I устанавливается равным некоторому постоянному значению ). Во-вторых, данное значение устанавливается не мгновенно, а в течение некоторого характерного времени т^). Не вдаваясь в детали динамики ионных каналов, для моделирования данных свойств можно использовать простое линейное дифференциальное уравнение следующего вида [13]:

(I ) - I

dt t(V)

(8)

Впервые модель, использующая такой упрощенный способ описания динамики ионных каналов нервных клеток, была введена в работах ФитцХью [14,15], а позже была обобщена Нагумо на случай пространственно-распределенных систем [16]. В рамках модели Фитцью-Нагумо используется подход, основанный на разделении различных ионных токов по их характерным временным масштабам (ср. с [17]). Предполагается, что характерные времена установления токов для большинства ионных каналов достаточно малы, что позволяет считать суммарный ток через них I^ всегда равным своему равновесному значению f (V). Выделяется также один или несколько типов каналов I, характерное время изменения токов которых maxy t(V) сравнимо и считается достаточно большим. Для описания этих токов используется уравнение типа (8), а для динамики трансмембранного потенциала записывается уравнение типа (5). В результате получим следующую обобщенную двумерную систему ФитцХью-Нагумо:

{ V = Iext - I - f (V), { I = e(g(V) - I),

где e = С/ maxt(V). Используя вид функций f и g из оригинальной работы Фитц-

Хью, после замены переменных получим классический вид системы ФитцХью-Нагумо

(и = и - и3/3 - v + Iext,

{ (10)

[ V = е(и - a - bv).

Здесь и описывает динамику мембранного потенциала нейрона; v - совокупное действие всех медленных ионных токов, отвечающих за восстановление потенциала покоя мембраны. Параметры a и b определяют проводимостные характеристики ионных каналов, а е (е > 0) - относительную скорость изменения медленных ионных токов.

Изначально на параметры b и е накладывались следующие условия: 0 < b < 1, е ^ 1/b. Первое условие обеспечивает существование в системе (10) лишь одного состояния равновесия, отвечающего состоянию покоя нейрона. Второе условие, с одной стороны, обеспечивает возможность смены этим состоянием равновесия устойчивости, а с другой - приводит к тому что, динамика системы в целом становится близкой к релаксационной. Таким образом, вне зависимости от параметров системы, ее динамика формируется так называемыми быстрыми и медленными движениями [18-21]. Первые происходят вдоль кривых v ~ const. Вторые происходят в е-окрестности устойчивого многообразия медленных движений, образованного двумя падающими участками и-нульклины (v = и - и3/3 + Iext). Левый падающий участок (в области отрицательного мембранного потенциала) обычно ассоциируют с рефрактерным (или невозбужденным) состоянием нейрона, а правый падающий участок (в области положительного мембранного потенциала) - с возбужденным состоянием нейрона.

Рассмотрим более подробно типичные режимы классической системы ФитцХью-Нагумо. Их качественные фазовые портреты приведены на рис. 5, a-f а на рис. 5, g-i приведены характерные отклики системы в этих режимах. Режим на рис. 5, a реализуется при a > a+. Состояние равновесия O в этом режиме является глобально устойчивым и расположено на левом участке устойчивого многообразия медленных движений. При этом действие стимула может вызывать генерацию системой одиночных импульсов возбуждения (потенциалов действия или спайков), поэтому этот режим называют еще возбудимым. Действительно, при достаточно сильной стимуляции изображающая точка (траектория Со) может покинуть окрестность состояния равновесия и преодолеть пороговое многообразие, которое образовано слоем неустойчивых медленных движений, происходящих в е-окрестности возрастающего участка и-нульклины, и его следом в обратном времени (на рис. 5, a область отмечена серым цветом). Далее по быстрым движениям изображающая точка переходит на правый участок устойчивого многообразия медленных движений, движется вдоль него до точки срыва, а затем по быстрым движениям вновь возвращается на левый участок, релаксируя в конце концов обратно к состоянию покоя. При a = a+ происходит бифуркация Андронова-Хопфа и смена устойчивости состояния равновесия, которое переходит на неустойчивое многообразие медленных движений. При уменьшении a < a+ на фазовой плоскости нейрона появляется малый предельный цикл Ci (см. рис. 5, b), отвечающий квазигармоническим подпороговым колебаниям. Далее система проходит через каскад так называемых уточных решений [20], содержащих как участки многообразий быстрых и устойчивых медленных движений, так

Рис. 5. Типичные режимы классической системы ФитцХью-Нагумо: (a, e, f) - режимы одиночной и периодической генерации импульсов возбуждения (потенциалов действия или спайков); (bd) - режимы подпороговых колебаний; (g-i) - характерные отклики в соответсвующих режимах. Пунктиром изображены нульклины системы. Параметры: (а) - a > a+; (b-f) - a £ (a-,a+), aT (s, b, Iext) = Ul,2 - b(ui,2 - uf,2/3 + Iext), Ui,2 =

Fig. 5. Typical regimes of classical FitzHugh-Nagumo system: (a, e, f) - regimes of single and periodic génération of excitation pulses (action potentials or spikes); (b-d) - regimes of subthreshold oscillations; (g-i) - typical responses in the regimes. The dashed lines represent the system's nullclines. The parameters:

(a) - a> a+; (b-f) - a £ (a-,a+), aT(e,b, Iext) = ui,2 - b(ui,2 - и^/З + Iext), ui,2 = TV 1 - bs

и участки неустойчивых многообразий медленных движений. Эти решения отвечают сложным подпороговым (см. рис. 5, c, d) и спайковым (см. рис. 5, e) колебаниям системы. Сначала малый предельный цикл деформируется в так называемую «малую утку без головы» (С2). Называют это решение так, потому, что оно содержит лишь один левый участок устойчивого многообразия медленных движений. «Малая утка» затем деформируется в «большую утку без головы» (C3), а та, в свою очередь, в «большую утку с головой» (C4). Последняя содержит как левый, так и правый участки устойчивого многообразия медленных движений. Наконец, «большая утка с головой» деформируется в большой предельный цикл C5 (см. рис. 5, f). Нейрон при этом переходит в режим периодической генерации импульсов возбуждения. При дальнейшем уменьшении a наблюдается обратный каскад «уточных» решений, при a = a- происходит бифуркация Андронова-Хопфа, состояние равновесия снова становится устойчивым, переходит на правый участок устойчивого многообразия медленных движений, и в системе возникает возбудимый режим.

Классическая система ФитцХью-Нагумо широко применяется и в настоящее время, причем в самых разнообразных областях современной науки. Однако ряд свойств этой системы затрудняют ее применение для описания динамики нейронов.

Пороговое многообразие формируется е-слоем неустойчивых медленных движений в окрестности неустойчивого медленного многообразия, и достаточно трудно однозначно классифицировать факт возбуждения. Основной режим динамики нейронов, возбудимый режим, реализуется в системе при достаточно малых е. Для преодоления этих недостатков систему рассматривают в области, где существует три состояния равновесия и пороговое многообразие формируется сепаратрисами седлового состояния (см. ниже), либо же разрабатывают на ее основе модифицированные системы. Модификации в основном касаются вида функции д в уравнении для медленных ионных токов. Например, в [22,23] теоретически и экспериментально изучалась модель с нелинейным восстановлением, отличающаяся кусочно-линейным видом этой функции,{

(и = и - и3/3 - V + 1ех1, (11)

) V = е(аиН(-и) + виН(и) - V - а),

где Н - функция Хевисайда, а параметры а, в (а > 0, в > 0) и а описывают нелинейные свойства ионных каналов. С одной стороны, такой вид функции д, более адекватно описывает динамику ионных токов, делая ее схожей с динамикой токов в системе Моррис-Лекара. С другой стороны, такая модификация приводит к возникновению новых динамических свойств системы. В частности, исследования показали, что в такой системе наблюдается большее по сравнению с классической моделью разнообразие динамических режимов.

Система (11) в зависимости от параметров может иметь либо одно, либо три грубых состояния равновесия. В первом случае это может быть как устойчивый, так и неустойчивый фокус (или узел). И в этой области параметров система ведет себя аналогично поведению классической системы ФитцХью-Нагумо, уже рассмотренному выше. Во втором случае на фазовой плоскости имеется два состояния равновесия типа фокус (или узел), которые могут быть как устойчивыми, так и неустойчивыми, и одно состояние равновесия седлового типа. На рис. 6 представлены разбиения плоскости параметров (е, 10 = а + 1ехг) для трех различных соотношений параметров а и в на области, отвечающие различным динамическим режимам системы (11). Фазовые портреты режимов, существующих соответственно в областях 1-15, приведены на рис. 7. Фазовые портреты в областях со штрихами зеркально эквивалентны соответствующим портретам в областях без штрихов. Кривые А+ и А+А-) на диаграммах отвечают, соответственно, бифуркациям Андронова-Хопфа (смены устойчивости) состояний равновесия О1 и Оз, локализованных в левой (и < 0) и правой (и > 0) полуплоскости (см. рис. 7). Верхний индекс при А1 и А3 показывает знак первой ляпуновской величины. Кривые Н+ и Н+(Н—) соответствуют бифуркациям «малых» петель сепаратрис седлового состояния равновесия - и - Щи, а кривые Н+2(Н-2) и Н+1 - бифуркациям «больших» петель сепаратрис и -Щи. Верхний индекс при Н1} Н2, Н12 и Н21 показывает знак седло-вой величины. Кривая С отвечает бифуркации образования «большого» двукратного предельного цикла, охватывающего все три состояния равновесия системы (11), а кривые С3а и С3 соответствуют бифуркациям образования «малых» двукратных предельных циклов, локализованных в окрестности состояния равновесия О3. Заметим, что разбиение на рис. 6 а можно получить и для классической системы ФитцХью-Нагумо [24], однако разбиения на рис. 6 Ь, с можно получить уже только в рамках моделей типа Моррис-Лекара.

0.8 ■

0.7

0.6

0.5

0.4

\ч "2+ 7 J

\ з\

2 N. 5 // ^Ъ^Ш/2'

С >Ч9 /7+

// 9 Н2,\ у/ 4 1

/Я!

0.8

0.6

0.4-

0.2-

---—_ 7

1,//: '"

-Щ,2 Ч

2 Щш ki

с..........■■¡J/J

Я^Ш

!/1

рлх

-0.078 0 0.022 0.044 /0 0 0.05 0.10 0.15 0.20 /п

а

7

сь гз

"г /Г 14

+ С Яй 8

Н\,2 >^5

i

0.885 0.880 0.875 0.870 0.865 0.860

с 0.0544 0.0560 0.0576

Рис. 6. Бифуркационные диаграммы различных динамических режимов системы (11), построенные для: a - а = 0.8, в = 0.9; b -а = 0.5, в = 2.0; с - а = 0.8, в = 1.1258

Fig. 6. Bifurcation diagrams for system (11) obtained for: a - а = 0.8, в = 0.9; b - а = = 0.5, в = 2.0; с - а = 0.8, в = 1.1258

В зависимости от параметров система может иметь до трех различных аттракторов. Так по одному аттрактору (моностабильность) имеется в областях 1 и 2, по два аттрактора (бистабильность) - в областях 3-5, 7-10 и по три аттрактора (мульти-стабильность) - во всех остальных областях. К аттракторам относятся: устойчивые состояния равновесия 0\ и Оз, отвечающие состояниям покоя и «перманентного» возбуждения нейрона, соответственно; глобальный устойчивый предельный цикл -Ья, охватывающий все три состояния равновесия и отвечающий режиму периодической генерации импульсов возбуждения; локальные устойчивые предельные циклы - ¡1 и ¡1* (¡| и ¡|*), охватывающие состояние равновесия О1 (03) и отвечающие режимам подпороговых колебаний (перманентного нестационарного возбуждения). Бассейны притяжения этих аттракторов ограничиваются устойчивыми сепаратрисами седла О2 и/или глобальным неустойчивым предельным циклом Ьи, и/или локальными неустойчивыми предельными циклами Iи, Щ, охватывающими состояния равновесия 01 и Оз, соответственно. Наличие тех или иных аттракторов и пороговых многообразий, а также их определенных конфигураций свидетельствует о наличии в системе определенного типа поведения. Если на фазовой плоскости имеются устойчивые предельные циклы, то имеет место колебательное поведение. Ес-

т п о

Рис. 7. Структурно устойчивые фазовые портреты системы (11), отвечающие динамическим режимам в областях 1-15 (см. рис. 6), соответственно

Fig. 7. Structurally stable phase portraits of system (11) corresponding to regimes in regions 1-15 (see fig. 6), respectively

ли обе неустойчивые сепаратрисы седла стремятся к одному и тому же состоянию равновесия, то имеет место возбудимое поведение относительно этого состояния. Если неустойчивые сепаратрисы стремятся к различным аттракторам или на фазовой плоскости имеются неустойчивые предельные циклы, то имеет место триггер-ное (переключательное) поведение. В зависимости от начальных условий в системе реализуется один из аттракторов, который впоследствии может меняться посредством соответствующей внешней стимуляции, перебрасывающей систему в область притяжения необходимого аттрактора. Заметим, что все рассмотренные выше типы поведения могут наблюдаться одновременно. Так, в области 1 наблюдается лишь колебательное поведение. В областях 2-6 наблюдается возбудимое поведение. Кроме того, в областях 3, 6 возможно установление или переход в состояние перманентного возбуждения, а в областях 3-6 возможен переход в различные колебательные режимы. Наконец в областях 7-15 наблюдается всевозможное триггерное поведение между состоянием покоя и перманентного возбуждения (области 7, 11-13), а также между различными колебательными режимами (области 8-15).

Одной из основных черт режимов, описанных выше (даже достаточно сложных), является наличие лишь не более чем двух порогов возбуждения, отделяющих друг от друга различное динамическое поведение системы (колебания, покой, возбудимость). Рассмотрим ещё ряд режимов системы (11), обладающих достаточно нетривиальными сложно-пороговыми динамическими свойствами, а именно, наличием значительного числа порогов возбуждения.

На рис. 8, a приведен фазовый портрет одного из таких режимов, который имеет место в области 1, при приближении параметров системы либо к бифуркационной кривой, отвечающей «большой» петле сепаратрис с отрицательной седловой величиной Н-2, либо к кривой «большого» двукратного предельного цикла С (Рис. 6). При этом входящая сепаратриса седла приближается к выходящей сепаратрисе ШЦ и начинает совершать вращение (при Ь ^ —ж) вокруг всех трех состояний равновесия. На фазовой плоскости (и, V) это соответствует появлению «слоений», ограниченных витками сепаратрисы . Число таких «слоений» растет с приближением параметров к обозначенным бифуркационным кривым и устремляется в бесконечность при точном совпадении параметров с бифуркационными значениями. При попадании начальных условий в одно из «слоений», возвращаясь в состояние равновесия О 1, траектория будет совершать некоторое число оборотов, что означает генерацию серии импульсов возбуждения. На рис. 8, Ь представлены области р¿, определяющие в зависимости от параметра е диапазон начальных возмущений ир (по переменной и при фиксированном V), отвечающих генерации различного числа импульсов (г = 0, 7). Область ро определяет отсутствие отклика, область р^ отклик в виде серии из г последовательных импульсов. Следует отметить, что генерация импульсов в данном режиме становится возможной не только для возбуждающих (ир > 0), но и подавляющих (ир < 0) воздействий. Рис. 8, ^ d иллюстрируют отклик системы на возбуждающее и подавляющее воздействия в виде 2 и 6 импульсов, соответственно. Подобный режим может возникать и в случае когда система имеет лишь одно состояния равновесия. Например, он наблюдается, когда параметры системы приближаются к бифуркационным значениям, отвечающим рождению в её фазовом пространстве двукратного предельного цикла. А образование «слоений» связано со сгущением фазовых траекторий, являющимся следствием динамической памяти системы о недавнем существовании двукратного предельного цикла.

Рис. 8. a - возбудимый мультипороговый режим. b - области генерации серий импульсов отклика (pi) в зависимости от величины внешнего воздействия и параметра е. c и d - генерация серий импульсов отклика в ответ на возбуждающее и подавляющее воздействия. Параметры: а = 0.8, в = 0.9,10 = 0.01

Fig. .8. a - excitable multi-threshold regime. b - areas of generation (pi) of different number of response pulses (spike trains) depending on the amplitude of external stimulus and the parameter е. c и d - the examples of responses to excitatory and inhibitory stimula. The parameters are а = 0.8, в = 0.9,10 = 0.01

Другой мультипороговый режим имеет место в области 7 в окрестности бифуркационной кривой «большого» двукратного предельного цикла C (см. рис. 6, a). Типичный фазовый портрет этого режима приведен на рис. 9, a. Здесь, как и в предыдущем случае, при приближении параметров к C, начинают вращаться обе входящие сепаратрисы седла. При этом триггерный бистабильный режим (переходы между состояниями равновесия) приобретает возбудимые мультипороговые свойства. Появляется, с одной стороны, возможность генерации импульсов и серий импульсов возбуждения причем как относительно нижнего, так и относительно верхнего состояния равновесия. С другой стороны, вращение устойчивых (входящих) сепаратрис приводит к возникновению нетривиального триггерного поведения, переходы между состояниями равновесия которыв также сопровождаются генерацией серий импульсов возбуждения. На рис. 9, b представлены области, отвечающие разнообразному поведению системы (11), возникающему в ответ на возмущение состояния покоя O\. Здесь области Si отвечают триггерному поведению системы (11), сопровождающемуся генерацией серии из i импульсов возбуждения. Рис. 9, c, d иллюстрируют типичное поведение системы (11) в данном режиме.

Еще один мультипороговый режим имеет место в области 10 также в окрестности бифуркационной кривой «большой» петли сепаратрис с отрицательной седло-вой величиной H-2 (см. рис. 6, c). Типичный фазовый портрет этого режима при-

Рис. 9. a - возбудимо-триггерный мультипороговый режим. b - области генерации серий импульсов отклика с возвращением к исходному (p¿) или переходом в другое (s¿) состояние равновесия. c и d -примеры откликов. Параметры: а = 0.8, в = 0.9, /о = 0.035

Fig. 9. a - excitable-trigger multi-threshold regime. b - areas of generation of different number of response pulses (spike trains) with either return to previous equilibrium state (pi) or transition to another one (s¿). c и d - the examples of responses. The parameters are а = 0.8, в = 0.9,10 = 0.035

веден на рис. 10, а. Здесь, в отличие от предыдущего случая, при приближении параметров к И-2, вокруг всех состояний равновесия начинают вращаться обе входящие сепаратрисы седла W( и W|. Это приводит к тому, что ранее простой режим с триггерно-подобной динамикой (переходы между состоянием покоя и локальными колебаниями) приобретает возбудимо-колебательное мультипороговое поведение. Пусть изначально (11) система находится в состоянии покоя. Тогда в зависимости от его начального возмущения можно получить генерацию серий импульсов возбуждения, заканчивающуюся либо возвращением обратно в состояние покоя, либо переходом в режим низкоамплитудных колебаний. Аналогичные процессы происходят при возмущении колебательного состояния системы. На рис. 10, b представлены области, отвечающие различному отклику системы (11) на возмущение ее состояния покоя. Как и прежде, область pi отвечает отклику системы в виде серии из i импульсов возбуждения с возвращением обратно в состояние покоя. Область же Ci также соответствует генерации серии из i импульсов, но с последующим переходом в колебательный режим. Рис. 10, c, d иллюстрируют оба рассмотренных типа поведения.

Аналогичными мультипороговыми свойствами обладает система (11) в области 8 для параметров, локализованных вблизи кривой C, и в областях 13 и 15 для параметров, локализованных вблизи кривой Н—. Фазовые портреты данных режимов представлены на рис. 11. Данные режимы имеют общие черты с рассмотренными

Рис. 10. a - возбудимо-колебательный мультипороговый режим. b - области генерации серий импульсов отклика (pi) с возвращением к исходному состоянию равновесия либо (ct) переходом в режим периодической генерации импульсов. c и d - примеры откликов. Параметры: а = 0.8, в = 1.1258,10 = 0.0595

Fig. 10. a - excitable-oscillatory multi-threshold regime. b - areas of generation of different number of response pulses (spike trains) with either return to previous equilibrium state (pt) or transition to regime of periodic spiking (ci). c и d - the examples of responses. The parameters are а = 0.8, в = 1.1258,10 = 0.0595

Рис. 11. (a, b) - возбудимые мультипороговые режимы с фоновым состоянием равновесия (O3) и предельным циклом (¿3), соответственно. c - возбудимо-колебательный мультипороговый режим с фоновым предельным циклом (13)

Fig. 11. (a, b) - excitable multi-threshold regimes with background equilibrium state (O3) and limit cycle (¿3), respectively. c - excitable-oscillatory multi-threshold regimes with background limit cycle (¿3)

выше режимами с мультипороговыми свойствами, поэтому останавливаться на них подробно не будем.

Подводя итог, можно сказать, что модель ФитцХью-Нагумо, представляющая собой простую динамическую систему на плоскости, адекватно описывает ключевые черты динамики нейрона - наличие потенциала покоя, возможность генерации потенциала действия, порог возбуждения, рефрактерность и проч. Наблюдаемые в данной системе динамические режимы включают в себя возбудимость и периодическую генерацию импульсов действия. Однако, в отличие от более сложных для аналитического исследования моделей типа Ходжкина-Хаксли, динамика модели ФитцХью-Нагумо может быть изучена качественными методами и проиллюстрирована на фазовой плоскости. Кроме того, важной особенностью моделей ФитцХью-Нагумо является достаточно простая схема радиотехнической реализации [16,23,24]. Именно благодаря простоте, с одной стороны, и адекватному отражению базовых свойств нейрона, с другой, модель ФитцХью-Нагумо пользуется большой популярностью среди исследователей.

2.4. Модель Хиндмарша-Роуза. Приведённые выше модели нейронной активности описывают генерацию нейронами отдельных импульсов - потенциалов действия, или спайков. Однако в нейрофизиологии известна другая, не менее важная, форма нейронной активности - генерация бёрстов. Бёрсты представляют собой группы двух или более спайков, идущих подряд друг за другом и перемежаемых периодами отсутствия активности. Впервые модель нейронного бёрстинга была предложена в работе Хиндмарша и Роуза [25]. Данная модель может трактоваться, с одной стороны, как обобщение модели ФитцХью-Нагумо, а с другой стороны как упрощение физиологически реалистичной модели Ходжкина-Хаксли. Данная модель является оптимальной с точки зрения баланса простоты и богатства демонстрируемых динамических режимов. Модель Хиндмарша-Роуза представляет собой трехмерную динамическую систему следующего вида:

y — ax3 + bx2 + Iext - z, (12)

c — dx2 — y, (13)

r (s (x — xi) — z). (14)

Здесь переменная x описывает мембранный потенциал нейрона, переменные y и z соответствуют быстрому и медленному ионным токам, протекающим через мембрану, причём параметр r ^ 1 определяет отношение характерных временных масштабов изменения этих токов. Параметр Iext описывает внешний ток, подаваемый на нейрон, а остальные параметры характеризуют нелинейность проводимости мембраны. Типичными значениями являются a = c = 1, b = 3, d = 5, остальные параметры обычно рассматриваются как контрольные и позволяют управлять динамическими режимами системы.

Благодаря малости параметра r в системе Хиндмарша-Роуза могут быть выделены так называемые быстрые и медленные движения [26]. Первые два уравнения для переменных x и y задают быструю подсистему, переменная z является медленной. При r = 0 значение переменной z является фиксированным и входит в быструю подсистему в качестве параметра. Наиболее важным свойством быстрой подсисте-

dx dt dy dt dz

dt

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

При 0 < г ^ 1 происходит медленное изменение переменной г. При этом аттракторы быстрой системы формируют в фазовом пространстве полной системы два многообразия медленных движений: состоянию равновесия в быстрой системе соответствует одномерное многообразие Ш0, а предельному циклу - двумерное многообразие представляющее собой осцилляционную поверхность (рис. 12, Ь). Устойчивые многообразия медленных движений системы по одной разделены неустойчивым многообразием Ши, соответствующим седлу быстрой подсистемы. Граница многообразия определяется гомоклинической бифуркацией в быстрой подсистеме, а граница многообразия Ш0 определяется седло-узловой бифуркацией в быстрой подсистеме.

В соответствии с уравнением (14), динамика медленной переменной г зависит от состояния быстрой подсистемы: если быстрая система находится в колебательном режиме, то есть на многообазии то переменная г растет, а если быстрая система находится в состоянии покоя, то есть на многообразии Шо, то медленная переменная г убывает. Таким образом в полной системе возникает повторяющееся движение по многообразиям. Сначала изображающая точка движется по многообразию в сторону увеличения г, пока не достигает его края. Затем система переключается на

Рис. 12. a - характерный фазовый портрет быстрой подсистемы (12)-(13) модели Хиндмарша-Роуза. В широкой области параметров система бистабильна, в ней сосуществуют устойчивое состояние равновесия O и устойчивый предельный цикл L, разделённые седлом U. b - качественное изображение полного фазового пространства модели Хиндмарша-Роуза. Устойчивые многообразия W0 и W1 разделены неустойчивым многообразием Wu. Также изображены три различных положения плоскости нульклины для уравнения (14)

Fig. 12. a - a typical phase portrait of the fast subsystem (12)-(13) of the Hindmarsh-Rose model. The system is bistable in a wide range of parameters with a stable steady state O and a stable limit cycle L separated by a saddle U. b - a qualitative view of the full phase space of the Hindmarsh-Rose model. The stable manifolds W0 and Wi are separated by the unstable manifold Wu. Also depicted are three different positions of the nulcline plane for Eq. (14)

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

Заметим, что скорость движения системы вдоль многообразий медленных движений определяется расположением плоскости z = s(x — xi), которая определяет нульклину уравнения (14). Положение данной плоскости контролируется параметрами s и xi, которые таким образом определяют длительности бёрстов и периодов затишья между ними. Так, при больших xi, когда плоскость нульклины пересекает многообразие Wi, система не покидает данного многообразия (положение (1) на рис. 12, b). В этом случае периодов затишья не возникает, и система генерирует периодические спайки с высокой частотой. При малых xi, когда плоскость нульклины пересекает седловое многообразие, происходят повторяющиеся переключения меду многообразиями (положение (3) на рис. 12, b). При этом чем ниже расположена плоскость нульклины, тем больше времени система проводит на многообразии Wo и, соответственно, тем дольше периоды затишья между бёрстами.

Большой интерес представляет промежуточное положение плоскости нулькли-ны, когда пересечение нульклины с многообразием Wi происходит вблизи седла (положение (2) на рис. 12, b). В данной области параметров происходит переход от генерации спайков к генерации бёрстов, ставший предметом ряда исследований [27-29]. Следует отметить, что в данной области параметров разделение движений системы на быстрые и медленные становится некорректным. Действительно, период предельного цикла в быстрой подсистеме неограниченно возрастает по мере приближения к точке бифуркации. Исследования показывают, что при уменьшении xi предельный цикл, соответствующий периодическому спайкингу, претерпевает каскад бифурка-

Рис. 13. Характерная динамика полной системы Хиндмарша-Роуза. а - последовательные переходы между многообразиями W0 и W1. b - зависимость мембранного потенциала x от времени - периодическая генерация бёрстов, состоящих из четырех спайков

Fig. 13. Typical dynamics of the Hindmarsh-Rose system. a - sequential transitions between the manifolds Wo and Wi. b - the membrane potential x versus time - periodic generation of burst consisting of fours spikes

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

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

2.5. Модель нейрона с постдеполяризацией. Как было сказано выше, важнейшим свойством нервных клеток является возбудимость, то есть способность генерировать электрический импульс в ответ на внешнее воздействие достаточной амплитуды. Количественной характеристикой возбудимости является порог возбуждения нейрона, то есть минимальная величина амплитуды внешнего воздействия, приводящее к возбуждению системы. Оказывается, порог возбуждения некоторых нейронов не является постоянной величиной, а изменяется во времени: генерация одного потенциала действия может привести к временному изменению порога возбуждения последующих импульсов. Данное явление называется постгиперполяриза-цией (айегЬурегро1ап2а1:юп, АНР) в случае повышения порога и постдеполяризацией (айегёеро1аг12а1:юп, АБР) в случае его понижения (см., например, [30-32]).

Свойство постдеполяризации нейронной мембраны может играть важную роль в динамике отдельных нейронов и нейронных ансамблей. Понижение порога возбуждения после однократной генерации потенциала действия может приводить к возникновению в дальнейшем спонтанных вторичных спайков или к переходу нейрона в режим бёрстовой активности. Предполагается, что постдеполяризация играет важную роль в когнитивных процессах центральной нервной системы, в частности, при формировании рабочей памяти [33-35]. Ионные механизмы, лежащие в основе постдеполяризации, до конца не изучены, однако наиболее вероятный механизм связан с некоторым ионным током, который быстро активируется после возбуждения нейрона и затем медленно инактивируется [36,37].

Свойство постдеполяризации описывает модель, предложенная в работе [38]. Данная система является обобщением классической модели ФитцХью-Нагумо, в которую вводится дополнительный входящий ионный ток г. Этот ток быстро включается после возбуждения нейрона и медленно уменьшается после прекращения этого возбуждения. Модель имеет следующий вид:

¿х

= г + 1ех* - У - ^(х), (15)

¿¿У = а + д(х) - Ьу, (16)

¿7

- = а(х)(1 - г) - в(х)г. (17)

Здесь х описывает мембранный потенциал, 1ехь - внешний ток, у и г соответствуют двум различным типам ионных токов, малый параметр е ^ 1 соответствует емкости мембраны, f (х) является нелинейной функцией кубической формы, д(х) -

монотонно возрастающая функция, функции а(х) и |3(х) описывают активацию и инактивацию входящего ионного тока z и имеют вид ступенчатых функций.

В системе (15)-(17) также можно выделить быстрые и медленные движения. Многообразие медленных движений системы изображено на рис. 14. Оно имеет размерность 2 и состоит из двух устойчивых и одной неустойчивой компонент. Устойчивая компонента Wo располагается в области отрицательных х и соответствует невозбужденному состоянию нейрона, компонента Wi расположена при положительных х и соответствует возбуждению.

Как и классическая модель ФитцХью-Нагумо, система (15)-(17) является возбудимой: под действием достаточно сильного импульса изображающая точка переходит с многообразия Wo на многообразие Wi, нейрон генерирует спайк. После этого изображающая точка возвращается на многообразие Wo и движется вдоль него. Важно, что в процессе движения траектория приближается к границе данного многообразия. В случае, если траектория не достигает этой границы, нейрон возвращается в состояние покоя (рис. 14, a). Однако важное значение имеет тот факт, что пока траектория проходит вблизи границы, мембранный потенциал х испытывает положительное отклонение по сравнению с равновесным состоянием. Это значит, что мембрана имеет в это время пониженный порог возбудимости, то есть для повторного возбуждения достаточно стимула меньшей амплитуды. Именно в этом и заключается свойство постдеполяризации (рис. 15, a). Динамический механизм, лежащий в основе данного свойства, заключается в приближении траектории, соответствующей потенциалу действия, к граничному множеству, определяющему порог возбуждения системы.

В случае, если траектория пересекает границу многообразия W0, нейрон повторно генерирует спайк, как это изображено на рис. 14, b. Этот процесс может повториться несколько раз или вылиться в периодическое возбуждение нейрона. Возникновение режима периодической генерации потенциалов действия связано с бифуркацией двойного предельного цикла. В результате этой бифуркации в системе рождается устойчивый предельный цикл, и система становится би-стабильной. Внешний стимул приводит к «переключению» нейрона из состояния покоя в режим периодической генерации потенциалов действия (рис. 15, b). Такое поведение может быть интерпретировано как динамическая память: нейрон «запо-

\ Z"

а

Рис. 14. Движения между многообразиями медленных движений: a - одиночный спайк; b - повторяющиеся спайки. Модифицировано из Klinshov V.V., Nekorkin V.I. // Communications inNonlinear Science and Numerical Simulation. 2012. Vol. 17(3). Pp. 1438-1446

Fig. 14. Transitions between the slow motion manifolds: a - a single spike; b - tonic spiking. Modified from Klinshov V.V., Nekorkin V.I. Communications in Nonlinear Science and Numerical Simulation, 2012, 17(3), pp.1438-1446

Рис. 15. Отклики нейрона на внешние стимулы, параметры е = 10~2, b = 0.5, So = 0.1, Si = 1, a = —1 зафиксированы, параметр у изменяется: a - возбудимость и постдеполяризация, у = 1.6; b - периодическая генерация спайков, у = 2. Модифицировано из Klinshov V.V., Nekorkin V.l. // Communications in Nonlinear Science and Numerical Simulation. 2012. Vol. 17(3). Pp. 1438-1446

Fig. 15. The response of a neuron on external stimuli, parameters е = 10~2, b = 0.5, So = 0.1, S i = 1, a = —1 are fixed, parameter у is varied: a - excitability and afterdepolarization, у = 1.6; b - tonic spiking, Y = 2. Modified from Klinshov V.V., Nekorkin V.I. Communications in Nonlinear Science and Numerical Simulation, 2012, 17(3), pp. 1438-1446

минает» внешнее воздействие и хранит информацию о нём в форме периодической активности (ср. с [33]).

Заметим, что вблизи точки бифуркации в системе могут наблюдаться так называемые спонтанные вторичные спайки [39,40]. Под действием внешнего импульса изображающая точка попадает в область сгущения траекторий и совершает несколько оборотов, после чего вновь возвращается в состояние равновесия. Это соответствует генерации нейроном нескольких дополнительных спайков после прекращения стимуляции.

2.6. Модель нейрона нижних олив. Нижние оливы являются одним из основных элементов системы контроля и координации движений [41] у многих позвоночных животных и наиболее ярко выражены в организмах птиц и млекопитающих. Нижние оливы специфическим образом преобразуют сигналы, поступающие от коры больших полушарий и спинного мозга, и передают эти сигналы далее в различные отделы мозжечка. Сигналы, поступающие через нижние оливы несут информацию о планируемом к выполнению движении, на основе которой мозжечок оптимизирует управление мышцами, обеспечивая рациональность, эффективность и отточенность движений.

Индивидуальная динамика нейронов нижних олив характеризуется ритмической колебательной активностью. Эта активность может проявляться как в форме квазипериодических подпороговых колебаний, так и в регулярной генерации потенциалов действия (спайков) [42,43]. Кроме того, нейроны нижних олив могут также демонстрировать хаотическую динамику.

Электрофизиологические свойства нейронов нижних олив определяются ионными токами, протекающими через клеточную мембрану [44]. Эти токи формируются в основном ионами кальция, калия и натрия и обладают различными уровнями активации и временами действия. Кальциевые ионные токи участвуют преимущественно в медленных процессах деполяризации и гиперполяризации клеточной мембраны и отвечают за генерацию низкочастотных подпороговых колебаний мембранного потенциала. Натриевые ионные токи отвечают за формирование коротких потенциалов действия, возникающих, как правило, на волне плавной деполяризации и/или в ответ на возбуждающее внешнее воздействие. Однако спайковый отклик может также воз-

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

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

du „, ч ^ dt =f (u) -

^ = U - (z - Ica) - iNa,

f

dz d(kt)

dw k d(kt)

(18)

= f(z) - w,

= eca{z - Ica - Iext(t)) ,

где переменные z и w отвечают за генерацию подпороговых колебаний и низкопороговых (Са2+-зависимых) спайков, а переменные u и v генерацию высокопороговых (Na2+-зависимых) спайков. Параметры еоа и £Na контролируют временные масштабы собственных колебаний модели, токи Ica и управляют уровнем деполяризации мембраны двух подсистем, f - нелинейная функция кубической формы (f (x)=x(x-a)(1-x)). Параметр к устанавливает отношение временных масштабов подсистем модели. Величина Iext(t) моделирует действие внешних токов на нейрон.

Заметим, что связь кальциевой и натриевой подсистем в (18) является односторонней. Иными словами изменение величины кальциевых ионных токов влечет за собой изменение натриевых, но не наоборот. Подобная конфигурация носит название «master-slave configuration». Таким образом, в модели (18) колебания кальциевых ионных токов задают ритм, на фоне которого формируется динамика натриевых ионных токов.

При соответствующем выборе параметров модель (18) может генерировать устойчивые подпороговые или спайковые колебания. Рассмотрим более подробно механизм возникновения таких колебаний, выбрав в качестве контрольного параметра величину Ic а. Положим вначале, что параметр Ica достаточно мал и при Iext(t) = 0 система (18) находится в состоянии покоя. С ростом Ica в кальциевой подсистеме происходит рождение периодических колебаний в результате суперкритической бифуркации Андронова-Хопфа. Эти колебания, в свою очередь, порождают вынужденные квазигармонические колебания в подсистеме для натриевых ионных токов. Такое поведение и определяет возникновение в модели нейрона нижних олив устойчивых подпороговых колебаний мембранного потенциала. Результаты численного моделирования демонстрируют хорошее качественное совпадение с данными натурных экспериментов (рис. 16), полученных in vitro при послойных исследовани-

Рис. 16. Сравнение результатов моделирования с экспериментальными данными, полученными in vitro. Значения параметров: eNa = 0.001,

еса = 0.02, iNa = —0.6, Ica = 0.01, a = 0.01,

к = 10. Модифицировано из [47]

Fig. 16. Comparison of simulation results with experimental data obtained in vitro. Parameter values: eNa = 0.001, eCa = 0.02, INa = —0.6, ICa = 0.01, a = 0.01, к = 10. Modified from [47]

Рис. 17. Сравнение результатов моделирования с экспериментальными данными, полученными in vitro. Значения параметров: eNa = 0.001, eCa = 0.02, INa = —0.59, ICa = 0.03, a = 0.01, к = 10. Модифицировано из [47]

Fig. 17. Comparison of simulation results with experimental data obtained in vitro. Parameter values: eNa = 0.001, eCa = 0.02, INa = = —0.59, ICa = 0.03, a = 0.01, к = 10. Modified from [47]

ях мозгового ствола морских свинок и опубликованных впервые в работе ЬНпаз & Уагош [47].

Амплитуда подпороговых колебаний в (18) при увеличении бифуркационного параметра 1Са нарастает вначале плавно, и дальнейшее поведение модели зависит от величины Если параметр достаточно велик, то с ростом 1Са в натриевой подсистеме амплитуда колебаний мембранного потенциала превышает порог возбудимости, и модель начинает генерировать периодические релаксационные колебания большой амплитуды - так называемые Ма2+-зависимые спайки (рис. 17). Заметим, что возникновение спайков происходит на волне гиперполяризации мембраны нейрона, вызванной колебаниями кальциевых ионных токов, то есть полностью определяется в этом случае фазой подпороговых колебаний. Если величина достаточно мала, то амплитуда колебаний в кальциевой подсистеме недостаточна для генерации Ма2+-зависимых спайков. Однако кальциевая подсистема обладает собственным порогом возбудимости. При превышении этого порога с ростом 1са форма колебаний в кальциевой подсистеме трансформируется нелинейным образом и при последующем увеличении 1са практически не изменяется. Мембранный потенциал при этом совершает колебания, форма которых существенно отличается от квазигармонической формы подпороговых колебаний, а амплитуда остается значительно ниже амплитуды Ма2+-зависимых спайков. Такое поведение соответствует режиму генерации Са2+ -зависимых спайков, наблюдаемых в натурных экспериментах (рис. 18). Отметим также, что в зависимости от параметров модель (18) может качественно воспроизводить спонтанные или вынужденные колебания, наблюдаемые в нейронах нижних олив.

Модель (18) представляет значительный интерес для моделирование коллективной динамики ансамблей нейронов нижних олив. Объединенные в единую сеть

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

2.7. Модели типа «накопление-сброс». Перейдем к описанию моделей нейронов, не претендующих на детальное описание ионного транспорта через мембрану. Такие модели нацелены преимущественно на воспроизведение базовых динамических режимов и свойств нейронов без привязки к определяющим их физическим механизмам. С одной стороны, эти модели часто используются при исследовании больших нейронных ансамблей, так как, в силу своей низкой размерности, позволяют сэкономить вычислительные ресурсы. С другой стороны, простота подобных моделей часто позволяет получать аналитические результаты для крупномасштабных сетей, например, путём сведения их динамики к редуцированным среднеполевым уравнениям [48-57].

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

Рис. 18. Сравнение результатов моделирования с экспериментальными данными, полученными in vitro. Значения параметров: eNa = 0.001, £Ca = 0.02, iNa = -0.6, Ica = 0.03, a = 0.01, к = 10. Модифицировано из [47]

Fig. 18. Comparison of simulation results with experimental data obtained in vitro. Parameter values: eNa = 0.001, eCa = = 0.02, INa = -0.6, Ica = 0.03, a = 0.01, к = 10. Modified from [47]

C<dt = IeXt - g(Vm - V)'

(19)

где Уо - потенциал покоя, 1ехъ - подаваемый на нейрон внешний ток, д - проводимость утечки мембраны. Когда мембранный потенциал достигает некоторого порогового значения У^ > Уо, считается, что нейрон генерирует спайк, после которого его мембранный потенциал мгновенно устанавливается равным некоторому значению

УеэеЛ < У0.

Динамика такой системы проиллюстрирована на рис. 19, a. Под действием внешнего тока мембранный потенциал постепенно увеличивается и при достаточно интенсивном воздействии в некоторый момент достигает порогового значения. После этого генерируется спайк, и мембранный потенциал мгновенно уменьшается до минимального значения и снова постепенно нарастает. Благодаря характерному для подобных систем чередованию периодов длительного накопления и кратковременного сброса потенциала, модели такого типа часто называют «накопление-сброс».

Рис. 19. a - динамика системы (19) - зависимость от времени мембранного потенциала при заданном внешнем воздействии; спайк соответствует мгновенному изменению величины мембранного потенциала с Vth на Vreset. b - эквивалентная электрическая схема

Fig. 19. a - dynamics of the system (19) - time dependence of the membrane potential for a given external action; spike corresponds to an instantaneous change in the value of the membrane potential from Vth to Vreset. b - equivalent electrical circuit

Накопление соответствует подпороговой динамике мембранного потенциала, а сброс -генерации потенциала действия.

Модели подобного типа адекватно описывают динамику мембранного потенциала ниже порога возбуждения, воспроизводя её интегрирующие свойства. С физической точки зрения подпороговая динамика данной модели аналогична динамике линейного шунтированного конденсатора (рис. 19, b). В отсутствие внешнего сигнала мембранный потенциал равен потенциалу покоя V0. Поступление на вход ненулевого внешнего сигнала приводит к зарядке конденсатора, однако наличие конечного сопротивления вызывает ток утечки. Таким образом, короткий импульс тока вида Iext = Ab(t) приведет к экспоненциально затухающему отклику u(t) = u0+(A/C) exp (-tg/C). Вследствие линейности подпороговой динамики, подача на нейрон серии коротких импульсов вызовет отклик, являющийся линейной суперпозицией откликов на каждый из импульсов. Таким образом, модель (19) эффективно интегрирует входные сигналы на характерных интервалах времени т ~ g/C.

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

В зависимости от задачи в модели «накопление-сброс» могут вводиться различные дополнения для более адекватного описания тех или иных свойств нейрона. Так, для учёта свойства рефрактерности зачастую вводится следующее уточнение. После генерации спайка мембранный потенциал в течение некоторого фиксированного времени Tref удерживается равным Vreset, при этом нейрон не реагирует на внешние воздействия (см., например, [59]). Для учёта свойства адаптации частоты

проводимость мембраны д может полагаться изменяющейся во времени и описываться отдельным дифференциальным уравнением, например,

йд

Та Ть = "д'

причём при каждой генерации нейроном спайка проводимость д мгновенно увеличивается на некоторую величину Ад. Такая модификация приводит к возрастанию проводимости мембраны при генерации спайка и, как следствие, к постепенному снижению частоты спайков при фиксированном внешнем воздействии. Возможны и другие методы учёта адаптивных свойств нейрона, например, путём введения дополнительного медленного ионного тока, а также учёта нелинейности мембраны [60-62].

Большой популярностью пользуется так называемая квадратичная модель типа «накопление-сброс» [63-65]. В данной модели динамика мембранного потенциала описывается следующим квадратичным уравнением:

йУ

— = 1еХг + а(У - Уо)(У - Уш). (20)

При достижении мембранным потенциалом величины У^ нейрон генерирует спайк, и его потенциал мгновенно приобретает значение Угезег. Квадратичная модель обладает рядом дополнительных свойств по сравнению с линейной моделью. Например, при определённых значениях параметров в ней возможна бистабильность между режимом покоя и периодической генерацией спайков. Модель (20) является нормальной формой для так называемой возбудимости первого класса, связанной с седло-узловой бифуркацией [66].

3. Дискретные модели нейронной активности

В последнее время для исследования кооперативных явлений в крупномасштабных нейронных сетях возрастающий интерес вызывает новый класс систем с дискретным временем: системы взаимосвязанных точечных отображений. Модели в форме точечных отображений обладают рядом преимуществ по сравнению с моделями в форме дифференциальных уравнений. Например, если для воспроизведения колебательных свойств в системах дифференциальных уравнений требуется, как минимум, два измерения, а хаотического поведения - три, то в дискретном времени и тот, и другой тип динамики можно описать в рамках даже одномерного отображения. Это преимущество проявляется при моделировании сложных режимов активности даже отдельных нейронов и в значительной степени крупномасштабных нейронных цепей, состоящих из взаимодействующих между собой различных структурных единиц. Например, для моделирования режима хаотических спайк-бёрстовых колебаний (одного из важнейших режимов нейронной активности) с помощью систем с непрерывным временем потребуется, как минимум, трёхмерная нелинейная система дифференциальных уравнений. С другой стороны, известны [67-69] дискретные модели в форме двумерных точечных отображений, которые достаточно хорошо воспроизводят режим спайк-бёрстовых колебаний. Точечные отображения позволяют моделировать и многие другие режимы нейронной активности. Например, модель

Д. Киалво [70] позволяет имитировать режимы так называемой нормальной и супернормальной возбудимости и др. Модель Н. Рулькова имеет несколько модификаций [71-73], одна из которых настроена на воспроизведение различных спайковых и бёрстовых колебательных режимов, а другая способна генерировать так называемые подпороговые колебания, то есть колебания малой амплитуды ниже порога возбуждения. Модель Курбажа-Некоркина [74, 75] является достаточно универсальной и описывает многие режимы нейронной активности.

3.1. Модель Киалво. Модель Киалво в форме точечного отображения имеет следующий вид [70]:

Хп+1 = х2пеУп-Хп + I, (21)

Уп+1 = ауп - Ьхп + с. (22)

Здесь переменная х воспроизводит динамику мембранного потенциала; у - восстанавливающая переменная; а, Ь, с, I являются параметрами. Эти параметры управляют динамикой системы следующим образом: а - константа времени восстановления (а < 1); Ь определяет степень зависимости восстанавливающих процессов от уровня активности; с - постоянное смещение; параметр I характеризует действие ионных токов, инжектированных в нейрон. В зависимости от параметра I динамика модели может быть разнообразной. Ключевую роль в установлении того или иного режима играет взаимное расположение изоклин у = 1п Х—1 + х и х = -(1-^у+с. Первая из этих изоклин представляет N-образную кривую с максимумом и минимумом, которая определяет пороговые свойства модели. Вторая изоклина - прямая, положение которой управляет устойчивостью неподвижной точки отображения.

В работе [70] описаны различные моды нейронной активности, воспроизводимые системой. Остановимся подробнее на некоторых из них.

Нормальная возбудимость. В таком режиме после генерации потенциала действия (спайка) и периода рефрактерности, мембранный потенциал монотонно возрастает к значению, соответствующему состоянию покоя. В модели последнему соответствует единственный аттрактор - устойчивый узел. Если амплитуда стимула достаточна для преодоления порога, на фазовой плоскости формируется траектория, напрямую возвращающаяся в неподвижную точку. Эта траектория воспроизводит режим нормальной возбудимости.

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

Колебательная бистабильность. При I = 0.03 динамика отображения биста-бильна. На фазовой плоскости одновременно сосуществуют два аттрактора: устойчивый фокус и устойчивая замкнутая инвариантная кривая. Их области притяжения разделяются неустойчивой замкнутой инвариантной кривой. В зависимости от начальных условий, траектории могут притягиваться либо устойчивой точкой, либо устойчивой инвариантной кривой. Соответственно, в этом случае отображение воспроизводит либо режим супернормальной возбудимости, либо периодических спай-ковых колебаний.

Хаотические спайк-бёрстовые колебания. Показано, что при а = 0.89, Ь = 0.18, с = 0.28, I = 0.03 в системе существуют неустойчивая неподвижная точка и хаотический аттрактор в области 0+. Динамический механизм хаотических колебаний состоит в следующем. Поскольку параметр Ь достаточно мал, переменная у остается более или менее в «замороженном» состоянии. Если пренебречь слагаемым -Ьх в уравнении (22), то динамика переменной у будет задаваться одномерным отображением

Уп+1 = ауп + с (23)

Это отображение имеет устойчивую неподвижную точку у = с/(1 — а) = уо с относительно большим мультипликатором в = 0.89. Поэтому любая траектория отображения стремится к этой точке достаточно быстро. В этом случае динамика переменной х описывается одномерным отображением

Хп+1 = хП ехр(уо — х) + I. (24)

Численное изучение этого отображения свидетельствует о его хаотических свойствах. Таким образом, относительно медленное изменение у и относительно быстрое изменение х формируют хаотический аттрактор в системе Киалво. Этот аттрактор определяет режим хаотических спайк-бёрстовых колебаний.

1.0-0,5 0 0.5 1.0 1.5 2.0 2.5 f -1.0-0.5 0 0.5 1.0 1.5 2.0 2.5 у

с d

Рис. 20. a - динамика переменных модели Киалво для режима спайковой активности. b - отображение быстрой системы. c и d - фазовая плоскость с изображенными изоклинами. Взят из [68]

Fig. 20. a - waveforms for the Chialvo model in case of the spiking activity regime. b - the map of fast motions. c and d - phase plane with nullclines. From [68]

3.2. Модель Ижикевича. Модель Ижикевича изначально разработана в непрерывном времени, но её дискретизация с временным шагом 1 Ш8 даёт точечное отображение [76]

Уп+1 = Р [Уп, I - Пп] , (25)

,Пп + а ■ [ЬУп - Пп], Уп < 30, ип+1 = \ , , . чп (26)

упп + а, Уп > 30,

где

( |шт(0.04у2 + 6у + 140 +1, 30), у< 30, (2)

| с, у > 30.

Здесь у - мембранный потенциал в шУ, и - медленная восстанавливающая переменная. Параметр а является малым и разделяет временные масштабы. Параметр с определяет напряжение при переустановке. Медленная подсистема имеет особенности: когда генерируется спайк, медленная переменная испытывает резкий скачок, равный значению параметра а. Быстро-медленные свойства модели позволяют ей воспроизводить разнообразные типы колебательной активности, в частности, спайк-бёрстовые колебания (рис. 21, a). Заметим, что отображение быстрых движений есть ни что иное как модель «накопление-сброс». Это видно из представленного рис. 21, Ь, на котором изображена диаграмма Кёнигса-Ламерея для медленной переменной и. Динамика полной системы определяется двумерным отображением, поэтому исчерпывающую информацию о возможных динамических сценариях нам дает фазовая плоскость (рис. 21, с) с изображенными изоклинами. Изоклина, получаемая из быстрого отображения, является параболой

и = 0.04у2 + 5у + 140 + I, у < 30, (28)

которая смещается влево или вправо на фазовой плоскости, если значение I уменьшается или увеличивается, соответственно. Поскольку парабола имеет две ветви, то при значениях переменной и, лежащих правее вершины параболы, быстрая система имеет две неподвижные точки: устойчивую и неустойчивую. Соответственно, ветви изоклины на фазовой плоскости (и, у) также могут быть классифицированы как устойчивая (N5) и неустойчивая (^и). Вершина параболы соответствует точке седло-узловой бифуркации быстрой подсистемы. С другой стороны, изоклина, получаемая из медленной подсистемы, является наклонной прямой. Когда изображающая точка находится выше этой прямой, значение и возрастает, когда ниже, значение и убывает. Параметр Ь управляет наклоном прямой изоклины. Множества, отмеченные на рис. 21, с и й как соответствуют спайковым циклам быстрой подсистемы. Их разрывы определяют смены периодов спайков. Когда множество М^^ лежит целиком выше медленной изоклины, значение и увеличивается вдоль спай-ковых орбит до и = иь, где линия переустановки у = с пересекает Nu, формируя гомоклиническую орбиту быстрой подсистемы.

В модели Ижикевича может быть от нуля до двух неподвижных точек. Их устойчивость управляется значением внешнего тока I, как это показано на рис. 21, с, й. Разница между ними заключается в горизонтальном смещении Ns, Мфкея. Когда значение I достаточно мало, одна неподвижная точка глобально устойчива, а другая является седлом. Нейрон находится в возбудимом режиме.

Рис. 21. a - динамика переменных модели Ижикевича для режима бёрстовой активности. b - отображение быстрой системы. c и d - фазовая плоскость с изображенными изоклинами. Взят из [68]

Fig. 21. a - waveforms for the Izhikevich model in case of the bursting activity regime. b - the map of fast motions. c and d - phase plane with nullclines. From [68]

Возмущение, переводящее за неустойчивую ветвь, вызывает генерацию ответного потенциала действия. Когда значение I растет, многообразия смещаются вправо и теряют устойчивость через бифуркацию Неймарка-Сакера при

Insizh = 16.25 - 62.5b + 6.25 I b2 - ( - a)

(b2 - (—П

V (1 - a)V .

(1 - а)2

Это значение лежит в окрестности точки пересечения медленной изоклины и вершины Ns - . Интересно отметить, что 5 перестает быть глобально устойчивой неподвижной точкой при слегка меньшем значении I, поскольку бифуркация Неймарка-Сакера субкритическая. После бифуркации 5 становится неустойчивым фокусом ^2, и в системе теперь существуют две неустойчивые неподвижные точки. Траектория притягивается орбитой, соответствующей спайковым или бёрстовым колебаниям.

3.3. Модель Рулькова. Двумерное отображение, описывающее динамику модели Рулькова, записывается в виде [ 71]

ХП+1 - f (xn,yn) 1

Уп+1 - Уп - И&п + 1) + ИО

(29)

(30)

Хп+\ 2.0

1.0

0

-1.0

-2.0

/

/

J,/ Щри

Ps f 1 1 1 1 1

где хп - быстрая переменная, отражающая взаимодействие калиевых и натриевых «быстрых» ионных токов; уп - медленная переменная, моделирующая динамическое действие «медленных» ионных токов, например, кальциевого тока. Медленное изменение во времени переменной уп обусловливается малостью параметра р,, в большинстве случаев будем полагать ^ = 0.001. Как будет показано ниже, параметр о задает невозмущенное состояние модельного нейрона. Быстрое одномерное отображение модельной системы построено так, чтобы оно могло отвечать как за спайковую активность нейрона, так и за его равновесное состояние. Такое поведение может быть достигнуто благодаря использованию разрывной функции специального

вида /(х, у)

а/(1 — х) + у, х < 0, /(х,у) = \ а + у, 0 <х< а + у,

— 1, х > а + у.

(31)

Здесь а является контрольным параметром отображения. Зависимость /(х,у) от х, полученная при фиксированном значении у, приведена на рис. 22. Значения параметров а и у выбраны таким образом, чтобы проиллюстрировать возможность одновременного существования неподвижных точек хрз, хри и предельного цикла Р%, соответствующего спайковым колебаниям. Функция строится таким образом, что когда значения у возрастают или уменьшаются, график зависимости /(х, у) от х соответственно поднимается вверх или опускается вниз. Исключение составляет лишь третий интервал, заданный неравенством х > а + у, где значения функции /(х,у) остаются равными (—1).

Уравнение для медленной переменной уп задает нелинейную обратную связь и позволяет моделировать как некоторые свойства медленных переходных процессов, так и медленные колебания, формирующие пачки импульсов нейронной активности. При значениях а < 4.0 в зависимости от величины параметра о отображение либо находится в режиме устойчивого равновесия, либо генерирует непрерывную последовательность спайков. Частота генерации импульсов при этом возрастает с увеличением значений параметра о. При а > 4 динамика отображения позволяет генерировать пачки спайков. Этот режим наблюдается в интервале значений параметра о, который располагается между режимами непрерывной импульсной активности и стационарным состоянием. Режимы пачечной активности включают в себя как регулярные, так и хаотические режимы.

-3.0 -2.0 -1.0 0 1.0 2.0

Рис. 22. Структура функции быстрого одномерного отображения, полученная при а = 5.6 и фиксированном yn = y = —3.75. График функции f (x,y) изображен жирной сплошной линией. Пунктирная линия изображает суперустойчивый цикл Pk. Устойчивая и неустойчивая неподвижные точки одномерного отображения обозначены как xps и xpu, соответственно. Взят из [71]

Fig. 22. Structure of the fast 1D map for the Rulkov model at а = 5.6 and fixed yn = y = —3.75. The bold solid line corresponds to f (x,y). The dotted line depicts the superstable cycle Pk. The stable and unstable fixed points of the 1D map are marked by xps and xpu, respectively. From [71]

3.4. Модель Курбажа-Некоркина.

жений [74,75] следующего вида:

Рассмотрим систему точечных отобра-

Хп+1 = Xn + F(xn) - вн(xn - d) - Уп, Уп+1 = Уп + e(xn - J)■

(32)

Здесь переменная х качественно характеризует изменение мембранного потенциала клетки, у отвечает за совокупное действие ионных токов (так называемая, восстанавливающая переменная). Параметр е определяет скорость изменения переменной у, параметры в, d, ,1 контролируют форму генерируемого сигнала. Отметим, что модель основана на дискретной версии известной в нейродинамике системы ФитцХью-

Нагумо с кубической нелинейностью F(х) и дополнительно введенной ступенчатой функцией Хевисайда Н(х)

А x=J

w\ о jß/^ JN. %\!f w w

t=J„,

c=d

с =ш

X

F(x) = x(x - а)(1 - x), "l, x > 0, 0, x < 0.

H (x) =

(33)

Рис. 23. Качественный вид изоклин горизонтальных и вертикальных наклонов на фазовой плоскости отображения (32), (33). Взят из [77]

Fig. 23. Shown qualitatively are the nullclines of map (32), (33). From [77]

На рис. 23 качественно показаны изоклины горизонтальных (прямая х = 1) и вертикальных (кривая у = = F(х) — вН(х — d)) наклонов системы (32), (33) на фазовой плоскости (х,у).

3.4.1. Регулярные режимы активности. Одним из основных свойств нейронов, изначально находящихся в состоянии покоя, является их способность к генерации потенциала действия, при превышении некоторого порога, в результате действия внешнего стимула (свойство возбудимости). Состоянию покоя нейрона в системе (32), (33) отвечает устойчивая неподвижная точка О. Показано, что на фазовой плоскости при этом существуют два порога, фактически определяемых неустойчивыми инвариантными кривыми ШЦ и ШЦ (точнее, тонкими слоями, состоящими из медленных траекторий, локализованных в окрестности этих инвариантных кривых). При действии на систему стимула, которого недостаточно для преодоления первого порога возбуждения (ШЦ), генерации потенциала действия не происходит, образуется лишь ответный отклик малой амплитуды (рис. 24, a, Ь (1)). Если амплитуда стимула оказывается достаточной для преодоления второго порога (ШЦ), то траектория попадает в область притяжения устойчивой инвариантной кривой и описывает характерную кривую, оканчивающуюся в устойчивой неподвижной точке О. Такому поведению на фазовой плоскости отвечает режим генерации одиночного потенциала действия, или спайка (рис. 24, a, Ь (и)).

Другой важный режим, который наблюдается в системе (32), (33) - это подпо-роговые колебания. На фазовой плоскости ему соответствует устойчивая замкнутая инвариантная кривая С (рис. 24, с), рожденная при смене устойчивости неподвижной точки в результате бифуркации Неймарка-Сакера. Еще один регулярный нейронный режим, воспроизводимый в системе, - периодические спайковые колебания. Установлено, что одним из условий их появления является относительная малость параметра в. При этом становятся возможными перемещения между слоями медленных движений, локализованных в окрестности двух устойчивых кривых и

Рис. 24. Фазовые портреты (левый столбец) и соответствующие осциллограммы (правый столбец) регулярных режимов нейронной активности в модели (32), (33). Взят из [77]

Fig. 24. Phase portraits (left column) and corresponding waveforms (right column) of regular regimes of neural activity in model (32), (33). From [77]

, без изменения направления движения при прохождении прямой разрыва х = й. В результате на фазовой плоскости формируется устойчивая замкнутая инвариантная кривая С (рис. 24, e), соответствующая периодическим спайковым колебаниям (рис. 24, /).

В модели (32), (33) помимо указанных наблюдается иной тип спайковой активности. При этом на фазовой плоскости образуется разрывный аттрактор А (рис. 24, g), определяющий колебания данной формы (рис. 24, к).

3.4.2. Хаотические режимы активности. В системе (32), (33) могут существовать различные хаотические аттракторы. Один из таких аттракторов, А, представлен на рис. 25, а, а соответствующий режим нейронной активности на рис. 25, Ь. В этом режиме потенциалы действия генерируются, перемежаясь с подпороговыми колебаниями. Динамический механизм таких колебаний заключается в том, что в окрестности линии разрыва х = й инвариантная кривая Жи разделяет траектории системы (32), (33) на два потока. Первый поток состоит из траекторий, которые совершают движения в окрестности х = й. Второй поток образован траекториями, преодолевшими второй порог и движущимися в окрестности кривой . В результате такого разделения траектория хаотически переходит из одного потока в другой, формируя аттрактор А.

Один из наиболее важных динамических режимов, встречающихся в различных нейронных структурах, - это хаотические спайк-бёрстовые колебания. В си-

Рис. 25. Фазовые портреты (левый столбец) и соответствующие осциллограммы (правый столбец) хаотических режимов нейронной активности в модели (32), (33). Взят из [77]

Fig. 25. Phase portraits (left column) and corresponding waveforms (right column) of chaotic regimes of neural activity in model (32), (33). From [77]

стеме (32), (33) данный режим с различными характеристиками воспроизводится в широком диапазоне значений параметров. На фазовой плоскости ему отвечает релаксационный хаотический аттрактор A (рис. 25, с), соответствующие ему колебания приведены на рис. 25, d. В силу важности режима хаотических спайк-бёрстовых колебаний, подробнее остановимся на динамическом механизме их установления.

3.4.3. Спайк-бёрстовые колебания. Хаотические спайк-бёрстовые колебания условно можно разделить на две фазы - быструю и медленную (см. рис. 25, c, d) Рассмотрим, как происходит формирование таких двухмасштабных колебаний. Пусть сначала е = 0 и, следовательно, переменная у = y° = const играет роль параметра в уравнении (32). В зависимости от значения у0 отображение быстрых движений (32) может демонстрировать как регулярную, так и хаотическую динамику (рис. 26). В случае регулярной динамики единственным аттрактором отображения (32) является устойчивая неподвижная точка xs (рис. 26, a, b). При квазистатическом уменьшении параметра у0 состояние равновесия xs при некотором бифуркационном значении сливается с неустойчивой неподвижной точкой xu и исчезает (происходит касательная, или седло-узловая бифуркация).

В случае хаотической динамики отображение (32) действует подобно отображению Лоренца и имеет инвариантный интервал, содержащий хаотический аттрактор (рис. 26, c, d). При квазистатическом увеличении параметра у0 этот аттрактор претерпевает внутренние бифуркации и, наконец, при некотором значении у0 = ycr происходит его граничный кризис [77]: одна из его границ сливается с неустойчи-

Рис. 26. Динамика отображения (32), (33) при различных значениях y0: a - единственная устойчивая неподвижная точка; b - устойчивая и неустойчивая неподвижные точки; c - устойчивая и неустойчивая неподвижные точки и инвариантный интервал, содержащий хаотический аттрактор; d - инвариантный интервал, содержащий хаотический аттрактор. Взят из [77]

Fig.26. Dynamics of map (32), (33) for different values of y0: a - a single stable fixed point; b - stable and unstable fixed points; c - stable and unstable fixed points and an invariant interval containing a chaotic attractor; d - an invariant interval containing a chaotic attractor. From [77]

вой неподвижной точкой хи, и аттрактор разрушается. Пусть теперь е > 0. В этом случае семейство устойчивых неподвижных точек отображения (32) образует на фазовой плоскости (х, у) при х < ■] устойчивую инвариантную кривую Wf. А при х > ■] семейство одномерных хаотических аттракторов образует переходное хаотическое множество. Допустим, что в начальный момент времени траектория стартует в окрестности кривой WS. Переменная у медленно убывает вдоль этой кривой, а переменная х поддерживается в квазиравновесном состоянии. Данная фаза образует медленный регулярный участок двумерного хаотического аттрактора, соответствующий пассивной фазе спайк-бёрстовых колебаний.

В окрестности значения у, при котором устойчивая Wl и неустойчивая WU кривые сливаются, траектория покидает слой медленных движений и попадает в область притяжения переходного хаотического множества. Поскольку теперь мы находимся в той части плоскости, где х > ■], то в этом случае переменная у начинает медленно расти. Итерируя отображение (32) численно, мы обнаружили, что хаотические колебания продолжают существовать в системе (32), (33) и после значения у0 = Усг, соответствующего кризису хаотического аттрактора в статическом случае е = 0 (см. рис. 26). Следовательно, имеет место эффект задержки исчезновения хаотических колебаний. Данная фаза движений образует быстрый хаотический участок аттрактора, который соответствует активной фазе или фазе деполяризации спайк-бёрстовых колебаний.

Отметим, что указанная универсальность модели Курбажа-Некоркина позволяет использовать данную модель при моделировании разнообразных процессов в нейроподобных системах [78-91].

Заключение

В обзоре представлены наиболее популярные динамические модели нейронов. Исторически первыми были сформулированы модели с непрерывным временем, такие как модель Ходжкина-Хаксли. Эта модель детально описывает ионный транспорт через нейронную мембрану, но является достаточно сложной с математической точки зрения. Наряду с ней представлены более простые редуцированные модели типа Моррис-Лекара, ФитцХью-Нагумо, Хиндмарша-Роуза и «накопление-сброс», которые, однако, обладают значительно меньшей степенью биологической релевантности. Кроме того, приводятся модели нейронов со специфическими свойствами, такими как последеполяризация и двухпороговая генерация спайков на фоне под-пороговых колебаний. Отдельную нишу занимают модели нейронов с дискретным временем, то есть в форме точечных отображений. С одной стороны, подобные модели не могут быть выведены из первых принципов и предлагают лишь качественное описание динамики реальных нейронов. С другой стороны, дискретные модели обладают большим богатством динамических режимов и удобны для численного моделирования, в том числе крупномасштабных нейронных сетей.

Работа выполнена в рамках государственного задания ИПФ РАН (проект № 0035-2014-0007) и программы РАН «Нелинейная динамика в математических и физических науках» (№ 0035-2018-006) при поддержке РФФИ (гранты № 17-02-00874, № 17-02-00904 и № 18-02-00406).

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

1. McCulloch W., Pitts W. A logical calculus of the ideas immanent in nervous activity // Bulletin of Mathematical Biophysics. 1943. Vol. 5, № 4. P. 115.

2. Hodgkin A.L., Huxley A.F. A quantitative description of membrane current and its application to conduction and excitation in nerve // J. Physiol. 1952. Vol. 117, № 4. P. 500.

3. Noble В. A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials // J. Physiol. 1962. Vol. 160, № 2. P. 317.

4. Plant R.E., Kim M.Mathematical description of a bursting pacemaker neuron by a modification of the Hodgkin-Huxley equations // Biophys. J. 1976. Vol. 16, № 3. P. 227.

5. Braun H.A., Huber M.T., DewaldM., Schafer K., and VoigtK. Computer simulations of neuronal signal transduction: The role of nonlinear dynamics and noise // Int. J. Bifurcation Chaos in Appl. Sci. Eng. 1998. Vol. 8. P. 881.

6. Morris C., Lecar H. Voltage oscillations in the barnacle giant muscle fiber // Biophys. J. 1981. Vol. 35. P. 193.

7. Keynes R.D., Rojas E., Taylor R. E., Vergara /.Calcium and potassium systems of a giant barnacle muscle fibre under membrane potential control // J. Physiol. (Lond.) 1973. Vol. 229, № 2. P. 409.

8. Gutkin B.S., Ermetrout G.B. Dynamics of membrane excitability determine interspike interval variability: a link between spike generation mechanisms and cortical spike train statistics // Neural Computation. 1998. Vol. 10, № 5. P. 1047.

9. Rinzel J., Ermetrout G.B. Analysis of neural excitability and oscillations // Methods in Neuronal Modeling: From Ions to Networks (Eds. C. Koch, I. Segev). London: MIT Press. 1999. P. 251.

10. Ermetrout G.B., Terman D.H.Mathematical Foundations of Neuroscience. New York: Springer. 2010. 422 p.

11. Tsumoto K., Kitajima H., Yoshinaga T., Aihara K., Kawakami H. Bifurcations in Morris-Lecar neuron model // Neurocomputing. 2006. Vol. 69, № 4-6. P. 293.

12. Behdad R., BinczakS., Dmitrichev A.S., Nekorkin V.I., Bilbault J.M. Artificial electrical Morris-Lecar neuron // IEEE Trans. Neural Netw. Learn. Syst. 2015. Vol. 26, № 9. P. 1875.

13. Abbot L.F. A network of oscillators // J. Phys. A: Math. Gen. 1990. Vol. 23. P. 3835.

14. FitzHugh R. Thresholds and plateaus in the Hodgkin-Huxley nerve equations // J. Gen. Physiol. 1960. Vol. 43. P. 867.

15. FitzHugh R. Impulses and physiological states in theoretical models of nerve membranes // Biophysical Journal. 1961. Vol. 1. P. 445.

16. Nagumo J., Arimoto S., Yoshizawa S. An active pulse transmission line simulating nerve axon // Proc. IRE. 1962. Vol. 50. P. 2061.

17. Kepler T.B., Abbott L.F., Marder E. Reduction of conductance-based neuron models // Biol. Cybern. 1992. Vol. 66, № 5. P. 381.

18. Андронов А.А., Витт А.А., Хайкин С.Э. Теория колебаний // М.: Физматиз, 1959. 915 с.

19. Mischenko E.F., Kolesov Yu.S., Kolesov A.Yu, Rozov N.Kh. Asymptotic Methods in Singularly Perturbed Systems, Monographs in Contemporary Mathematics // NY, Consultants Bureau. 1984. 294 p.

20. Арнольд В.И., Афраймович В.С., Ильяшенко Ю.С., Шильников Л.П. Теория бифуркаций // Итоги науки и техники. Серия «Современные проблемы математики. Фундаментальные направления». 1986. Т. 5. 218 с.

21. Fenichel N. Geometric singular perturbation theory for ordinary differential equation // SIAM J. Diff. Eqns. 1979. Vol. 31. P. 53.

22. Некоркин В.И., Дмитричев А.С., Щапин Д.С., Казанцев В.Б. Динамика модели нейрона со сложно-пороговым возбуждением // Математическое моделирование. 2005. Т. 17. № 6. C. 75.

23. Binczak S., Kazantsev V.B., Nekorkin V.I., Bilbault J.M. Experimental study of bifurcations in a modified FitzHugh-Nagumo cell//Electron. Lett. 2003. Vol. 39. P. 13.

24. Щапин Д.С. Динамика двух нейроноподобных элементов с подавляющей обратной связью // Радиотехника и электроника. 2009. Т. 54. № 2. С. 185.

25. Hindmarsh J.L., Rose R.M. A model of neuronal bursting using three coupled first order differential equations // Proc. of the Royal Society London B. 1984. Vol. 221. P. 87.

26. Некоркин В.И.Лекции по основам теории колебаний. Н. Новгород: Изд-во НН-ГУ, 2012. 311 с.

27. Wang X-J.Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle // Physica D. 1993. Vol. 62, № 1-4. P. 263.

28. Innocenti G., MorelliA., Genesio R., TorciniA. Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos // Chaos. 2007. Vol. 17, № 4. P. 043128.

29. Shilnikov A., Kolomiets M. Methods of the qualitative theory for Hindmarsh-Rose model: A case study - A tutorial // Int. J. of Bifurcation and Chaos. 2008. Vol. 18, № 8. P. 2141.

30. Izhikevich E.M. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge: MIT Press. 2007. 441 p.

31. Miura R.M. Analysis of excitable cell models // Journal of Computational and Applied Mathematics. 2002. Vol. 144, № 1-2. P. 29.

32. Yue C., Remy S., Su H., Beck H., Yaari Y. Proximal persistent Na+ channels drive spike afterdepolarizations and associated bursting in adult CA1 pyramidal cells // J. Neurosci. 2005. Vol. 25, № 42. P. 9704.

33. Lisman J.E., Idiart M.A. Storage of 7 +/- 2 short-term memories in oscillatory subcycles // Science. 1995. Vol. 267, № 5203. P. 1512.

34. Jensen O., Idiart M.A.P. and Lisman J.E. Physiologically realistic formation of autoassociative memory in networks with theta/gamma oscillations: Role of fast NMDA channels // Learn. Mem. 1996. Vol. 3, № 2-3. P. 243.

35. Jensen O., Lisman J.E. Hippocampal sequence-encoding driven by a cortical multi-item working memory buffer // Trends in Neurosciences. 2005. Vol. 28, № 2. P. 67.

36. Haj-Dahmane S., Andrade К. Ionic mechanism of the slow afterdepolarization induced by muscarinic receptor activation in rat prefrontal cortex // J. Neurophysiol. 1998. Vol. 80, № 3. P. 1197.

37. Park J.-Y., Remy S., Varela О., Cooper D.C., Chung S., Kang H.-W., Lee J.-H., Spruston N. A post-burst after depolarization is mediated by group i metabotropic glutamate receptor-dependent upregulation of Ca(v)2.3 R-type calcium channels in CA1 pyramidal neurons // PLoS Biology. 2010. Vol. 8, № 11. P. e1000534.

38. Клиньшов В.В., Некоркин В.И. Модель нейрона с последеполяризацией и краткосрочная память // Изв. ВУЗов Радиофизика. 2005. Vol. 48, № 3. P. 228.

39. Kepler T.B., Marder E. Spike initiation and propagation on axons with slow inward currents // Biol. Cybern. 1993. Vol. 68, № 3. P. 209.

40. Enns-Ruttan J., Miura R.M.Spontaneous secondary spiking in excitable cells // J. Theor. Biol. 2000. Vol. 205, № 2. P. 181.

41. Schweighofer N., Lang E.J., Kawato M. Role of the olivo-cerebellar complex in motor learning and control // Front. Neural Circuits. 2013. Vol. 7. Art. № 94. P. 1.

42. Manor Y., Rinzel J., Segev I., Yarom Y. Low-amplitude oscillations in the inferior olive: A model based on electrical coupling of neurons with heterogeneous channel densities // J. Neurophysiol. 1997. Vol. 77, № 5. P. 2736.

43. Velarde M.G., Nekorkin V.I., Kazantsev V.B., Makarenko V.I., Llinas R. Modeling inferior olive neuron dynamics // Neural Netw. 2002. Vol. 15, № 1. P. 5.

44. Schweighofer N., Doya K., Kawato M. Electrophysiological properties of inferior olive neurons: A compartmental model // J. Neurophysiol. 1999. Vol. 82, № 2. P. 804.

45. Kazantsev V.B., Nekorkin V.I., Makarenko V.I., Llinas R. Olivo-cerebellar cluster-based universal control system // Proc. Natl. Acad. Sci. USA. 2003. Vol. 100, № 22. P. 13064.

46. Kazantsev V.B., Nekorkin V.I., Makarenko V.I., Llinas R. Self-referential phase reset based on inferior olive oscillator dynamics // Proc. Natl. Acad. Sci. USA. 2004. Vol. 101, № 52. P. 18183.

47. Llinas R., Yarom Y Oscillatory properties of guinea-pig inferior olivary neurones and their pharmacological modulation: An in vitro study // J. Physiol. 1986. Vol. 376. P. 163.

48. Klinshov V., Franovic I. Slow rate fluctuations in a network of noisy neurons with coupling delay // EPL (Europhysics Letters). 2016. Vol. 116, № 4. P. 48002.

49. Klinshov V., Franovic I. Mean field dynamics of a random neural network with noise // Physical Review E. 2015. Vol. 92, № 6. P. 62813.

50. Franovic I., Klinshov V. Clustering promotes switching dynamics in networks of noisy neurons // Chaos. 2018. Vol. 28. P. 23111.

51. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons // Journal of Computational Neuroscience. 2000. Vol. 8, № 3. P. 183.

52. Olmi S., Politi A., Torcini A. Collective chaos in pulse-coupled neural networks // EPL (Europhysics Letters). 2010. Vol. 92, № 6. P. 60007.

53. Ullner E., Politi A. Self-sustained irregular activity in an ensemble of neural oscillators // Physical Review X. 2016. Vol. 6, № 1. P. 1.

54. Hasegawa H. Population rate codes carried by mean, fluctuation and synchrony of neuronal firings // Physica A: Statistical Mechanics and its Applications. 2009. Vol. 388, №4. P. 499-513.

55. Hasegawa H. Synchrony and variability induced by spatially correlated additive and multiplicative noise in the coupled Langevin model // Physical Review E. 2008. Vol. 78, № 3. P. 31110.

56. Nykamp D.Q., Friedman D., Shaker S., Shinn M., Vella M., Compte A., Roxin A. Mean-field equations for neuronal networks with arbitrary degree distributions // Physical Review E. 2017. Vol. 95, № 4. P. 1.

57. Montbrio E., Pazo D., Roxin A. Macroscopic description for networks of spiking neurons // Physical Review X. 2015. Vol. 5, № 2. P. 1.

58. Lapicque L. Recherches quantitatives sur l'excitation electrique des nerfs traitee comme une polarisation // J. Physiol. Pathol. 1907. Gen. 9. P. 620.

59. LazarA.A. Time encoding with an integrate-and-fire neuron with a refractory period //Neurocomputing. 2004. Vol. 58. P. 53.

60. Liu Y.-H., Wang X-J.Spike-frequency adaptation of a generalized leaky integrate-and-fire model neuron // Journal of Computational Neuroscience. 2001. Vol. 10, № 1. P. 25.

61. Brette R., Gerstner W. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity // Journal of Neurophysiology. 2005. Vol. 94, № 5. P. 3637.

62. Abbot L.F., van Vreeswijk C. Asynchronous states in networks of pulse-coupled oscillators // Physical Review E. 1993. Vol. 48. P. 1483.

63. Ermentrout B. Type I membranes, phase resetting curves, and synchrony // Neural Computation. 1996. Vol. 8, № 5. P. 979.

64. Latham P.E., Richmond B.J., Nelson P., Nirenberg S. Intrinsic dynamics in neuronal networks. I. Theory // J. Neurophysiology. 2000. Vol. 83, № 2. P. 808.

65. Hansel D., Mato G. Existence and stability of persistent states in large neuronal networks // Phys. Rev. Letters. 2001. Vol. 86. P. 4175.

66. Izhikevich E.M. Neural excitability, spiking and bursting // International Journal of Bifurcation and Chaos. 2000. Vol. 10, № 6. P. 1171.

67. Courbage M., Nekorkin V.I. Map based models in neurodynamics // International Journal of Bifurcation and Chaos. 2010. Vol. 20, № 6. P. 1631.

68. Ibarz B., Casado J.M., Sanjuan M. A. F. Map-based models in neuronal dynamics // Physics Reports. 2011. Vol. 501, № 1-2. P. 1.

69. Girardi-Schappo M., Tragtenberg M.H.R., Kinouchi O. A brief history of excitable map-based neurons and neural networks // Journal of neuroscience methods. 2013. Vol. 220, № 2. P. 116.

70. Chialvo D.R. Generic excitable dynamics on a two-dimensional map // Chaos, Solitons & Fractals. 1995. Vol. , № 3-4. P. 461.

71. Rulkov N.F. Modeling of spiking-bursting neural behavior using two-dimensional map // Physical Review E. 2002. Vol. 65, № 4. P. 041922.

72. Shilnikov A.L., Rulkov N.F. Subthreshold oscillations in a map-based neuron model // Physics Letters A. 2004. Vol. 328, № 2-3. P. 177.

73. Rulkov N.F. Regularization of synchronized chaotic bursts // Physical Review Letters. 2001. Vol. 86, № 1. P. 183.

74. Некоркин В.И., Вдовин Л.В. Дискретная модель нейронной активности // Известия вузов. Прикладная нелинейная динамика. 2007. Том. 15, № 5. Стр. 36.

75. Courbage M., Nekorkin V.I., Vdovin L.V. Chaotic oscillations in a map-based model of neural activity // Chaos: An Interdisciplinary Journal of Nonlinear Science. 2007. Vol. 17, № 4. P. 043109.

76. Izhikevich E.M., Hoppensteadt F. Classification of bursting mappings // International Journal of Bifurcation and Chaos. 2004. Vol. 14, № 11. P. 3847.

77. Maslennikov O.V., Nekorkin V.I. // Nonlinear Dynamics and Complexity (Eds. V. Afraimovich, A.C.J. Luo, X. Fu). Springer, 2014. P. 143.

78. Hess A., Yu L., Klein I., De Mazancourt M., Jebrak G., Mal H., Brugiere O., Fournier M., Courbage M., Dauriat G. Neural mechanisms underlying breathing complexity // PloS one. 2013. Vol. 8, № 10. P. e75740.

79. Courbage M., Maslennikov O.V., Nekorkin V.I. Synchronization in time-discrete model of two electrically coupled spike-bursting neurons // Chaos, Solitons & Fractals. 2012. Vol. 45, № 5. P. 645.

80. Maslennikov O.V., Nekorkin V.I. Modular networks with delayed coupling: Synchronization and frequency control // Physical Review E. 2014. Vol. 90, № 1. P. 012901.

81. Maslennikov O.V., Nekorkin V.I. Discrete model of the olivo-cerebellar system: structure and dynamics // Radiophysics and Quantum Electronics. 2012. Vol. 55, № 3. P. 198.

82. Nekorkin V.I., Maslennikov O.V Spike-burst synchronization in an ensemble of electrically coupled discrete model neurons // Radiophysics and Quantum Electronics 2011. Vol. 54, № 1. P. 56.

83. Maslennikov O.V., Nekorkin V.I., Kurths /.Basin stability for burst synchronization in small-world networks of chaotic slow-fast oscillators // Physical Review E. 2015. Vol. 92, № 4. P. 042803.

84. Maslennikov O.V., Nekorkin V.I. Evolving dynamical networks with transient cluster activity // Communications in Nonlinear Science and Numerical Simulation. 2015. Vol. 23, № 1-3. P. 10.

85. Yu L., De Mazancourt M., Hess A., Ashadi F.R., Klein I., Mal H., Courbage M., Mangin L. Functional connectivity and information flow of the respiratory neural network in chronic obstructive pulmonary disease // Human Brain Mapping. 2016. Vol. 37, № 8. P. 2736.

86. Maslennikov O.V., Kasatkin D.V., Rulkov N.F., Nekorkin V.I. Emergence of antiphase bursting in two populations of randomly spiking elements // Physical Review E. 2013. Vol. 88, № 4. P. 042907.

87. Maslennikov O.V., Shchapin D.S., Nekorkin V.I. Transient sequences in a hypernetwork generated by an adaptive network of spiking neurons // Phil. Trans. R. Soc. A. 2017. Vol. 375, № 2096. P. 20160288.

88. Yue Y., Liu Y.-J., Song Y.-L., Chen Y., Yu L.-C. Information Capacity and Transmission in a Courbage-Nekorkin-Vdovin Map-Based Neuron Model // Chinese Physics Letters. 2017. Vol. 34, № 4. P. 048701.

89. Franovic I., Maslennikov O.V., Bacic I., Nekorkin V.I. Mean-field dynamics of a population of stochastic map neurons // Physical Review E. 2017. Vol. 96, № 1. P. 012226.

90. Mangin L., Courtage MRespiratory Neural Network: Activity and Connectivity // Advances in Dynamics, Patterns, Cognition (Eds. Aronson, I.S., Pikovsky, A., Rulkov, N.F., Tsimring, L.S.). Nonlinear Systems and Complexity. Springer Int. 2017. Vol. 20. P. 227.

91. Yang X., Wang M.The evolution to global burst synchronization in a modular neuronal network // Modern Physics Letters B. 2016. Vol. 30, № 14. P. 1650210.

References

1. McCulloch W., Pitts W. A logical calculus of the ideas immanent in nervous activity. Bulletin of Mathematical Biophysics, 1943, vol. 5, № 4, p. 115.

2. Hodgkin A.L., Huxley A.F. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 1952, vol. 117, № 4, p. 500.

3. Noble B. A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pacemaker potentials. J. Physiol., 1962, vol. 160, № 2, p. 317.

4. Plant R.E., Kim M. Mathematical description of a bursting pacemaker neuron by a modification of the Hodgkin-Huxley equations. Biophys. J., 1976, Vol. 16, № 3. p. 227.

5. Braun H.A., Huber M.T., Dewald M., Schafer K., and Voigt K. Computer simulations of neuronal signal transduction: The role of nonlinear dynamics and noise. Int. J. Bifurcation Chaos Appl. Sci. Eng., 1998, vol. 8, p. 881.

6. Morris C., Lecar H. Voltage oscillations in the barnacle giant muscle fiber. Biophys. J., 1981, vol. 35, p. 193.

7. Keynes R.D., Rojas E., Taylor R. E., Vergara J. Calcium and potassium systems of a giant barnacle muscle fibre under membrane potential control. J. Physiol. (Lond.), 1973, vol. 229, № 2, p. 409.

8. Gutkin B.S., Ermetrout G.B. Dynamics of membrane excitability determine interspike interval variability: a link between spike generation mechanisms and cortical spike train statistics. Neural Computation, 1998, vol. 10, № 5, p. 1047.

9. Rinzel J., Ermetrout G.B. Analysis of neural excitability and oscillations. Methods in Neuronal Modeling: From Ions to Networks (Eds. C. Koch, I. Segev). London: MIT Press. 1999. P. 251.

10. Ermetrout G.B., Terman D.H. Mathematical Foundations of Neuroscience. New York: Springer. 2010. 422 p.

11. Tsumoto K., Kitajima H., Yoshinaga T., Aihara K., Kawakami H. Bifurcations in Morris-Lecar neuron model. Neurocomputing, 2006, vol. 69, № 4-6, p. 293.

12. Behdad R., Binczak S., Dmitrichev A.S., Nekorkin V.I., Bilbault J.M. Artificial electrical Morris-Lecar neuron. IEEE Trans. Neural Netw. Learn. Syst., 2015, vol. 26, № 9, p. 1875.

13. Abbot L.F. A network of oscillators. J. Phys. A: Math. Gen., 1990, vol. 23, p. 3835.

14. FitzHugh R. Thresholds and plateaus in the Hodgkin-Huxley nerve equations. J. Gen. Physiol., 1960, vol. 43, p. 867.

15. FitzHugh R. Impulses and physiological states in theoretical models of nerve membranes. Biophysical Journal, 1961, vol. 1, p. 445.

16. Nagumo J., Arimoto S., Yoshizawa S. An active pulse transmission line simulating nerve axon. Proc. IRE., 1962, vol. 50, p. 2061.

17. Kepler T.B., Abbott L.F., Marder E. Reduction of conductance-based neuron models. Biol. Cybern., 1992, vol. 66, № 5, p. 381.

18. Andronov A.A., Vitt A.A., Khaikin S.E. Theory of Oscillators. Oxford: Pergamon Press, 1966. 815 p.

19. Mischenko E.F., Kolesov Yu.S., Kolesov A.Yu., Rozov N.Kh. Asymptotic Methods in Singularly Perturbed Systems, Monographs in Contemporary Mathematics. NY, Consultants Bureau, 1984. 294 p.

20. Arnold V.I., Afrajmovich V.S., Il'yashenko Yu.S., Shil'nikov L.P. Bifurcation theory and catastrophe theory. Encyclopaedia of Mathematical Sciences: Dynamical Systems V, 1994, vol. 5, 274 p.

21. Fenichel N. Geometric singular perturbation theory for ordinary differential equation. SIAM J. Diff. Eqns, 1979, Vol. 31, P. 53.

22. Nekorkin V.I., Dmitrichev A.S., Shapin D.S., Kazantsev V.B. Dynamics if a neuron model with complex-threshold excitation. Mathematical Models and Computer Sim., 2005, vol. 17, № 6, p. 75. (In Russian).

23. Binczak S., Kazantsev V.B., Nekorkin V.I., Bilbault J.M. Experimental study of bifurcations in a modified FitzHugh-Nagumo cell. Electron. Lett., 2003, vol. 39, p. 13.

24. Shchapin D.S. Dynamics of two neuronlike elements with inhibitory feedback. Journal of Communications Technology and Electronics, 2009, vol. 54, № 2, p. 175.

25. Hindmarsh J.L., Rose R.M. A model of neuronal bursting using three coupled first order differential equations. Proc. of the Royal Society London B., 1984, vol. 221, p. 87.

26. Nekorkin V.I. Introduction to Nonlinear Oscillations. Wiley-VCH, 2015. 264 p.

27. Wang X.-J. Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle. Physica D., 1993, vol. 62, № 1-4, p. 263.

28. Innocenti G., Morelli A., Genesio R., Torcini A. Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos. Chaos, 2007, vol. 17, № 4, p. 043128.

29. Shilnikov A., Kolomiets M. Methods of the qualitative theory for Hindmarsh-Rose model: A case study - A tutorial. Int. J. of Bifurcation and Chaos, 2008, vol. 18, № 8, p. 2141.

30. Izhikevich E.M. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge: MIT Press, 2007. 441 p.

31. Miura R.M. Analysis of excitable cell models. Journal of Computational and Applied Mathematics, 2002, vol. 144, № 1-2, p. 29.

32. Yue C., Remy S., Su H., Beck H., Yaari Y. Proximal persistent Na+ channels drive spike afterdepolarizations and associated bursting in adult CA1 pyramidal cells. J.Neurosci, 2005, vol. 25, № 42, p. 9704.

33. Lisman J.E., Idiart M.A. Storage of 7 +/-2 short-term memories in oscillatory subcycles. Science, 1995, vol. 267, № 5203, p. 1512.

34. Jensen O., Idiart M.A.P. and Lisman J.E. Physiologically realistic formation of

autoassociative memory in networks with theta/gamma oscillations: Role of fast NMDA channels. Learn. Mem., 1996, vol. 3, № 2-3, p. 243.

35. Jensen O., Lisman J.E. Hippocampal sequence-encoding driven by a cortical multi-item working memory buffer. Trends in Neurosciences, 2005, vol. 28, № 2, p. 67.

36. Haj-Dahmane S., Andrade K. Ionic mechanism of the slow afterdepolarization induced by muscarinic receptor activation in rat prefrontal cortex. J. Neurophysiol., 1998, vol. 80, №3, p. 1197.

37. Park J.-Y., Remy S., Varela O., Cooper D.C., Chung S., Kang H.-W., Lee J.-H., Spruston N. A post-burst after depolarization is mediated by group is metabotropic glutamate receptor-dependent upregulation of Ca(v)2.3 R-type calcium channels in CA1 pyramidal neurons. PLoS Biology, 2010, vol. 8, № 11, p. e1000534.

38. Klin'shov V.V., Nekorkin V.I. Model of a neuron with afterdepolarization and short-term memory. Radiophysics and Quantum Electronics, 2005, vol. 48, № 3, p. 203.

39. Kepler T.B., Marder E. Spike initiation and propagation on axons with slow inward currents. Biol. Cytern., 1993, vol. 68, № 3, p. 209.

40. Enns-Ruttan J., Miura R.M. Spontaneous secondary spiking in excitable cells. J. Theor. Biol., 2000, vol. 205, № 2, p. 181.

41. Schweighofer N., Lang E.J., Kawato M. Role of the olivo-cerebellar complex in motor learning and control. Front. Neural Circuits, 2013, vol. 7, art. № 94, p. 1.

42. Manor Y., Rinzel J., Segev I., Yarom Y. Low-amplitude oscillations in the inferior olive: A model based on electrical coupling of neurons with heterogeneous channel densities. J. Neurophysiol., 1997, vol. 77, № 5, p. 2736.

43. Velarde M.G., Nekorkin V.I., Kazantsev V.B., Makarenko V.I., Llinas R. Modeling inferior olive neuron dynamics. Neural Netw., 2002, vol. 15, № 1, p. 5.

44. Schweighofer N., Doya K., Kawato M. Electrophysiological properties of inferior olive neurons: A compartmental model. J. Neurophysiol., 1999, vol. 82, № 2, p. 804.

45. Kazantsev V.B., Nekorkin V.I., Makarenko V.I., Llinas R. Olivo-cerebellar cluster-based universal control system. Proc. Natl. Acad. Sci. USA, 2003, vol. 100, № 22, p. 13064.

46. Kazantsev V.B., Nekorkin V.I., Makarenko V.I., Llinas R. Self-referential phase reset based on inferior olive oscillator dynamics. Proc. Natl. Acad. Sci. USA, 2004, vol. 101, №52, p. 18183.

47. Llinas R., Yarom Y. Oscillatory properties of guinea-pig inferior olivary neurones and their pharmacological modulation: An in vitro study. J. Physiol., 1986, vol. 376, p. 163.

48. Klinshov V., Franovic I. Slow rate fluctuations in a network of noisy neurons with coupling delay. EPL (Europhysics Letters), 2016, vol. 116, № 4, p. 48002.

49. Klinshov V., Franovic I. Mean field dynamics of a random neural network with noise. Physical Review E., 2015, vol. 92, № 6, p. 62813.

50. Franovic I., Klinshov V. Clustering promotes switching dynamics in networks of noisy neurons. Chaos, 2018, vol. 28, p. 23111.

51. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. Journal of Computational Neuroscience, 2000, vol. 8, № 3, p. 183.

52. Olmi S., Politi A., Torcini A. Collective chaos in pulse-coupled neural networks. EPL (Europhysics Letters), 2010, vol. 92, № 6, p. 60007.

53. Ullner E., Politi A. Self-sustained irregular activity in an ensemble of neural oscillators. Physical Review X., 2016, vol. 6, № 1, p. 1.

54. Hasegawa H. Population rate codes carried by mean, fluctuation and synchrony of neuronal firings. Physica A: Statistical Mechanics and its Applications, 2009, vol. 388, № 4, p. 499.

55. Hasegawa H. Synchrony and variability induced by spatially correlated additive and multiplicative noise in the coupled Langevin model. Physical Review E., 2008, vol. 78, №3, p. 31110.

56. Nykamp D.Q., Friedman D., Shaker S., Shinn M., Vella M., Compte A., Roxin A. Mean-field equations for neuronal networks with arbitrary degree distributions. Physical Review E., 2017, vol. 95, № 4, p. 1.

57. Montbrio E., Pazo D., Roxin A. Macroscopic description for networks of spiking neurons. Physical Review X., 2015, vol. 5, № 2, p. 1.

58. Lapicque L. Recherches quantitatives sur l'excitation electrique des nerfs traitee comme une polarization. J. Physiol. Pathol. Generale, 1907, 9, p. 620.

59. Lazar A.A. Time encoding with an integrate-and-fire neuron with a refractory period. Neurocomputing, 2004, vol. 58, p. 53.

60. Liu Y.-H., Wang X.-J. Spike-frequency adaptation of a generalized leaky integrate-and-fire model neuron. Journal of Computational Neuroscience, 2001, vol. 10, № 1, p. 25.

61. Brette R., Gerstner W. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. Journal of neurophysiology, 2005, vol. 94, № 5, p. 3637.

62. Abbot L.F., van Vreeswijk C. Asynchronous states in networks of pulse-coupled oscillators. Physical Review E., 1993, vol. 48, p. 1483.

63. Ermentrout B. Type I membranes, phase resetting curves, and synchrony. Neural Computation, 1996, vol. 8, № 5, p. 979.

64. Latham P.E., Richmond B.J., Nelson P., Nirenberg S. Intrinsic dynamics in neuronal networks. I. Theory. J. Neurophysiology, 2000, vol. 83, № 2, p. 808.

65. Hansel D., Mato G. Existence and stability of persistent states in large neuronal networks. Phys. Rev. Letters, 2001, vol. 86, p. 4175.

66. Izhikevich E.M. Neural excitability, spiking and bursting. International Journal of Bifurcation and Chaos, 2000, vol. 10, № 6, p. 1171.

67. Courbage M., Nekorkin V.I. Map based models in neurodynamics. International Journal of Bifurcation and Chaos, 2010, vol. 20, № 6, p. 1631.

68. Ibarz B., Casado J.M., Sanjuan M.A.F. Map-based models in neuronal dynamics. Physics Reports, 2011, vol. 501, № 1-2, p. 1.

69. Girardi-Schappo M., Tragtenberg M.H.R., Kinouchi O. A brief history of excitable map-based neurons and neural networks. Journal of Neuroscience Methods, 2013, vol. 220, №2, p. 116.

70. Chialvo D.R. Generic excitable dynamics on a two-dimensional map. Chaos, Solitons & Fractals, 1995, vol. , № 3-4, p. 461.

71. Rulkov N.F. Modeling of spiking-bursting neural behavior using two-dimensional map. Physical Review E., 2002, vol. 65, № 4, p. 041922.

72. Shilnikov A.L., Rulkov N.F. Subthreshold oscillations in a map-based neuron model. Physics Letters A., 2004, vol. 328, № 2-3, p. 177.

73. Rulkov N.F. Regularization of synchronized chaotic bursts. Physical Review Letters, 2001, vol. 86, № 1, p. 183.

74. Nekorkin V.I., Vdovin L.V. Map-based model of the neural activity. Izvestiya VUZ. Applied Nonlinear Dynamics, 2007, vol. 15, № 5, p. 36. (In Russian).

75. Courbage M., Nekorkin V.I., Vdovin L.V. Chaotic oscillations in a map-based model of neural activity. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2007, vol. 17, №4, p. 043109.

76. Izhikevich E.M., Hoppensteadt F. Classification of bursting mappings. International Journal of Bifurcation and Chaos, 2004, vol. 14, № 11, p. 3847.

77. Maslennikov O.V., Nekorkin V.I. // Nonlinear Dynamics and Complexity (Eds. V. Afraimovich, A.C.J. Luo, X. Fu). Springer, 2014. P. 143.

78. Hess A., Yu L., Klein I., De Mazancourt M., Jebrak G., Mal H., Brugiere O., Fournier M., Courbage M., Dauriat G. Neural mechanisms underlying breathing complexity. PloS ONE, 2013, vol. 8, № 10, p. e75740.

79. Courbage M., Maslennikov O.V., Nekorkin V.I. Synchronization in time-discrete model of two electrically coupled spike-bursting neurons. Chaos, Solitons & Fractals, 2012, vol. 45, № 5, p. 645.

80. Maslennikov O.V., Nekorkin V.I. Modular networks with delayed coupling: Synchronization and frequency control. Physical Review E., 2014, vol. 90, № 1, p. 012901.

81. Maslennikov O.V., Nekorkin V.I. Discrete model of the olivo-cerebellar system: structure and dynamics. Radiophysics and Quantum Electronics, 2012, vol. 55, № 3, p. 198.

82. Nekorkin V.I., Maslennikov O.V. Spike-burst synchronization in an ensemble of electrically coupled discrete model neurons. Radiophysics and Quantum Electronics, 2011, vol. 54, № 1, p. 56.

83. Maslennikov O.V., Nekorkin V.I., Kurths J. Basin stability for burst synchronization in small-world networks of chaotic slow-fast oscillators. Physical Review E., 2015, vol. 92, № 4, p. 042803.

84. Maslennikov O.V., Nekorkin V.I. Evolving dynamical networks with transient cluster activity. Communications in Nonlinear Science and Numerical Simulation, 2015, vol. 23, № 1-3, p. 10.

85. Yu L., De Mazancourt M., Hess A., Ashadi F.R., Klein I., Mal H., Courbage M., Mangin L. Functional connectivity and information flow of the respiratory neural network in chronic obstructive pulmonary disease. Human Brain Mapping, 2016, vol. 37, № 8, p. 2736.

86. Maslennikov O.V., Kasatkin D.V., Rulkov N.F., Nekorkin V.I. Emergence of antiphase bursting in two populations of randomly spiking elements. Physical Review E., 2013, vol. 88, № 4, p. 042907.

87. Maslennikov O.V., Shchapin D.S., Nekorkin V.I. Transient sequences in a hypernetwork generated by an adaptive network of spiking neurons. Phil. Trans. R. Soc. A., 2017, vol. 375, № 2096, p. 20160288.

88. Yue Y., Liu Y.-J., Song Y.-L., Chen Y., Yu L.-C. Information capacity and transmission in a Courbage-Nekorkin-Vdovin map-based neuron model. Chinese Physics Letters, 2017, vol. 34, № 4, p. 048701.

89. FranoviC I., Maslennikov O.V., BaCiC I., Nekorkin V.I. Mean-field dynamics of a population of stochastic map neurons. Physical Review E., 2017, vol. 96, № 1, p. 012226.

90. Mangin L., Courbage M. Respiratory neural network: Activity and connectivity. Advances in dynamics, patterns, cognition (Eds. Aronson, I.S., Pikovsky, A., Rulkov, N.F., Tsimring, L.S.). Nonlinear Systems and Complexity. Springer Int. 2017. vol. 20. p. 227.

91. Yang X., Wang M. The evolution to global burst synchronization in a modular neuronal network. Modern Physics Letters B., 2016, vol. 30, № 14, p. 1650210.

Дмитричев Алексей Сергеевич - родился в 1980 году в селе Слободское Кстовского района Горьковской области. Окончил радиофизический факультет Нижегородского государственного университета им. Н.И. Лобачевского (2002). С 2002 года работает в Институте прикладной физики РАН. Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (2010, ННГУ). Область научных интересов - нелинейная динамика, динамические сети, структурообразование, волны и простанственно-временной хаос, синхронизация. Имеет более 50 научных публикаций по указанным направлениям.

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: admitry@neuron.appl.sci-nnov.ru

Касаткин Дмитрий Владимирович - родился в 1978 году в Горьком. Окончил радиофизический факультет Нижегородского государственного университета им. Н.И. Лобачевского (2001). С 2002 года работает в Институте прикладной физики РАН. Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2004, ННГУ). Область научных интересов -теория нелинейных колебаний и волн, теория синхронизации, динамический хаос, нейродинамика. Имеет более 50 научных публикаций по указанным направлениям.

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: kasatkin@appl.sci-nnov.ru

Клиньшов Владимир Викторович - родился в 1984 году в Горьком, окончил радиофизический факультет Нижегородского государственного университета им. Н.И. Лобачевского (2006). С 2003 года работает в Институте прикладной физики РАН. Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2009, ИПФ РАН). Область научных интересов - нелинейная динамика, динамические сети, синхронизация, системы с запаздыванием. Имеет более 60 научных публикаций по указанным направлениям. В 2017 году награжден Медалью Российской академии наук с премией для молодых ученых РАН.

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: vladimir.klinshov@gmail.com

Кириллов Сергей Юрьевич - родился в 1986 году в Горьком, окончил радиофизический факультет Нижегородского государственного университета им. Н.И. Лобачевского (2009). С 2008 года работает в Институте прикладной физики РАН. Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2017, ННГУ им. Н.И. Лобачевского). Область научных интересов - нелинейная динамика, нейродинамика, теория бифуркаций, параметрические и автоколебательные системы, нелинейные волны в активных средах. Имеет более 20 научных публикаций по указанным направлениям.

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: skirillov@ipfran.ru

Масленников Олег Владимирович - родился в 1988 году в Горьком, окончил радиофизический факультет Нижегородского государственного университета им. Н. И. Лобачевского (2011). С 2007 года работает в Институте прикладной физики РАН. Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2014, ИПФ РАН). Область научных интересов -нелинейные колебания и волны, синхронизация, сложные сети, нейродинамика. Имеет более 40 научных публикаций по указанным направлениям. Награжден медалью РАН с премией для молодых ученых (2016).

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: olmaov@ipfran.ru

Щапин Дмитрий Сергеевич - родился в 1981 году в Горьком, окончил радиофизический факультет Нижегородского государственного университета им. Н.И. Лобачевского (2003). С 2003 года работает в Институте прикладной физики РАН. Защитил диссертацию на соискание учёной степени кандидата физико-математических наук (2011, ННГУ). Область научных интересов - теория нелинейных колебаний и волн, синхронизация, нейродинамика. Имеет более 10 научных публикаций по указанным направлениям.

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: shapinds@mail.ru

Некоркин Владимир Исаакович - родился в 1948 году, окончил радиофизический факультет Горьковского государственного университета им. Н.И. Лобачевского (1971). Доктор физико-математических наук (1992), профессор кафедры теории колебаний и автоматического регулирования ННГУ, заведующий отделом нелинейной динамики Института прикладной физики РАН, лауреат премии им. А.А. Андронова (2012). Область научных интересов - динамика нелинейных систем, нейродинамика, теория синхронизации, пространственно-временной хаос, структурообразование и нелинейные волны. Имеет более 250 научных публикаций в отечественных и зарубежных изданиях, в том числе 4 книги и 2 патента.

Россия, 603950 Нижний Новгород, 603950, ул. Ульянова, 46 Институт прикладной физики РАН E-mail: vnekorkin@neuron.appl.sci-nnov.ru