Санкт-Петербург, г. Санкт-Петербург и Ленинградская область, Россия
с 01.01.2015 по 01.01.2019
Санкт-Петербург, Россия
Санкт-Петербург, г. Санкт-Петербург и Ленинградская область, Россия
Санкт-Петербург, г. Санкт-Петербург и Ленинградская область, Россия
ВАК 03.01.2002 Биофизика
УДК 62 Инженерное дело. Техника в целом. Транспорт
ГРНТИ 27.35 Математические модели естественных наук и технических наук. Уравнения математической физики
ГРНТИ 27.41 Вычислительная математика
Проведено исследование с целью повышения качества компьютерного моделирования нейроморфных систем за счет идентификации нежелательных численных эффектов. В качестве методов исследования использованы вычислительный эксперимент, компьютерный анализ и синтез алгоритмов разностных схем интегрирования. Описано влияние эффектов дискретизации на динамическое поведение биологических моделей нейронов, получены инструменты поиска таких эффектов.
модель нейрона, эффекты дискретизации, нелинейная система, численное интегрирование, детерминированный хаос, хаотический переходный процесс
Задачи по определению характера функциональных процессов в структурах головного мозга зачастую решаются путем исследования математических моделей нейронов и синапсов. Понимание режимов работы нейронных сетей необходимо при разработке инструментов анализа мозговой активности и выявления возможных патологий. К примеру, состояние, предшествующее эпилептическому припадку, связано с сильной синхронизацией нейронов в пораженной области [1]. Также известно, что шизофрения характеризуется аномальными нервными колебаниями и синхронизацией [2]. Таким образом, современный взгляд на природу импульсной активности нейронов направлен на теорию динамических систем. Ходжкин и Хаксли были первыми, кто изучил динамику нейронов в таком понимании, представив модель мембраны аксона нейрона в 1952 г. [3]. До сих пор модель ионных каналов Ходжкина - Хаксли (ХХ) является наиболее распространенным кинетическим представлением чувствительной к напряжению проводимости в нейронах. Практически все компактные математические модели, демонстрирующие хорошее соответствие экспериментальным данным, основаны на классической модели ХХ.
Компьютерное моделирование генерации и передачи нервных импульсов обычно проводится для крупных нейронных сетей. Этот процесс предполагает значительные вычислительные затраты в случае использования сложных феноменологических моделей нейронов и синапсов. По этой причине для экономии временного ресурса все еще применяются численные методы низких порядков точности. Кроме того, встречаются публикации, в которых результаты исследований в области нелинейных динамических систем приводятся без явного указания примененных численных методов. Такие случаи отрицательно сказываются на воспроизводимости результатов моделирования.
В основе компьютерного моделирования динамических систем лежит применение дискретных операторов, которые предположительно способны сохранить ключевые свойства непрерывного прототипа [4; 5]. При этом на дискретные модели негативно влияют эффекты, вызванные численным методом, шагом дискретизации и ошибками округления. Этот факт особенно важен при анализе нелинейных систем, в том числе моделей биологических нейронов. В рамках данного исследования рассматривается влияние распространенных одношаговых методов интегрирования путем проведения ряда вычислительных экспериментов, в которых модель нейрона, представленная классическими уравнениями ХХ, переводится в режимы резонансных и хаотических колебаний. Результаты исследования демонстрируются с помощью двухпараметрических динамических карт, межспайковых интервальных гистограмм и бифуркационных диаграмм.
Модель нейрона Ходжкина - Хаксли
Система уравнений Ходжкина - Хаксли [3] представляет собой классическую феноменологическую модель нейрона, устанавливающую связь между внутренними ионными токами и динамическим поведением электрически возбудимой мембраны. Применяя экспериментальный метод фиксации напряжения на гигантском аксоне кальмара, Ходжкин и Хаксли определили, что в клеточной мембране протекают три основных тока: открываемый напряжением устойчивый ток К+, открываемый напряжением устойчивый ток Na+ и постоянный ток утечки, переносимый ионами Cl-. Эквивалентная электрическая схема для классической модели ХХ показана на рис. 1. Эта схема включает емкость мембраны C, сопротивление утечки RL, две чувствительные к напряжению нелинейные проводимости GNa и GK, описывающие ионные каналы натрия и калия, а также три батареи ENa, EK и EL, соединенные последовательно резистивным элементам.
Рис. 1. Эквивалентная электрическая схема модели нейрона ХХ
Система уравнений ХХ для эквивалентной схемы имеет вид
(1)
где
В уравнениях системы (1) переменная n моделирует четыре активационных затвора (следовательно, множитель n4) для тока K+, переменная m представляет три активационных затвора (множитель m3), h является переменной одного инактивационного затвора для тока Na+.
Наблюдаемые численные эффекты
Компьютерное моделирование динамической системы (1) выполнено в среде NI LabVIEW 2018 на основе следующих методов численного решения обыкновенных дифференциальных уравнений (ОДУ) с фиксированным размером шага интегрирования h: явный метод Эйлера первого порядка (Explicit Euler, EE), явный метод средней точки второго порядка (Explicit Midpoint, EMP), явный метод Рунге-Кутты четвертого порядка (Runge-Kutta 4, RK4) и явный метод Дормана-Принса восьмого порядка (DOPRI 78). Результаты, полученные с использованием метода DOPRI 78, взяты за эталонное решение.
Модель нейрона ХХ может демонстрировать стационарное состояние, режимы затухающих подпороговых колебаний (рис. 2а) и генерации потенциала действия (рис. 2б). Следовательно, согласно классификации Ижикевича [6], модель ХХ является бистабильным резонатором.
Рис. 2. Колебательные режимы системы (1) во временной области: а - подпороговые колебания;
б - генерация потенциалов действия
Первая серия экспериментов была проведена для демонстрации различий в поведении дискретных операторов в ответ на входной ток I амплитуды A = 2 мкА/см2 с различными значениями периода T и длительности импульса τ. Параметры системы (1) следующие: C = 1 мкФ/см2, = 36 мСм/см2,
= 120 мСм/см2, GL = 0,3 мСм/см2, EK = -12 мВ, ENa = 115 мВ, EL = 10,6 мВ. Начальные условия при моделировании системы (1) - [0,31; 0,05; 0,59; 0].
На рис. 3 представлена бифуркационная диаграмма по параметру T при τ = 5,5 мс. В интервалах 11–12 мс и 16–18 мс дискретная модель, полученная явным методом Эйлера с шагом h = 0,05 мс, генерирует потенциалы действия, что является качественным отклонением от прочих моделей на диаграмме, находящихся в режиме подпороговых колебаний.
Рис. 3. Бифуркационная диаграмма для импульсного входного сигнала
Нагляднее данный эффект можно представить с помощью динамических карт на рис. 4, где области белого цвета соответствуют стационарному режиму и режиму подпороговых колебаний, а цветом выделены области режима генерации потенциалов действия для рассматриваемых дискретных операторов. Горизонтальной линией на рис. 4а отмечена строка карты, построенная по бифуркационной диаграмме на рис. 3.
Рис. 4. Динамические карты для импульсного входного сигнала: а - разница между методом
первого порядка и эталоном; б - разница между методом второго порядка и эталоном
Как видно из рис. 4а, явный метод Эйлера первого порядка существенно увеличивает область возбуждения по сравнению с эталонным решением, что может негативно сказаться при моделировании нейронных сетей даже в режимах гармонических колебаний. Поведение дискретной модели нейрона на базе метода второго порядка приближено к эталонному (рис. 4б).
Вторая серия вычислительных экспериментов была проведена для случая, когда динамическая система (1) переходит в хаотический режим генерации потенциалов действия в ответ на сигнал I = I0 + A sin(2πωt). Параметры системы (1) остались прежними, за исключением ENa = 120 мВ. Начальные условия установлены на [0,31; 0,05; 0,59; 0,001].
В качестве инструмента анализа были выбраны межспайковые интервальные гистограммы. Диаграммы этого типа отображают количество интервалов n между максимумами потенциалов действия относительно их продолжительности t и обычно применяются при изучении статистических особенностей нейронов в ответ на зашумленный входной сигнал. На рис. 5 представлен ряд гистограмм (150 интервалов) для дискретных моделей в ответ на входной ток I0 = 6,22 мкА с добавлением 70 Гц синусоидальной компоненты амплитудой A = 0,6 мкА. Время моделирования установлено от 1000 до 3000 мс.
Из рис. 5 видно, что модели на основе EMP, RK4 и DOPRI предоставляют качественно схожие решения, где разброс интервалов свидетельствует о хаотическом режиме на всем временном промежутке моделирования. Дискретная модель, полученная методом EE, на данном промежутке оказалась в периодическом режиме. Анализ данной модели во временной области выявил наличие хаотического переходного процесса, длительность которого чувствительно зависит от начальных условий. Это явление подробно описано в работе [7].
Рис. 5. Межспайковые интервальные гистограммы
Похожий эффект был также найден для модели, полученной методом четвертого порядка RK4 при параметрах входного сигнала I0 = 6,48 мкА, A = 0,56 мкА и ω = 70 Гц. На рис. 6 показаны результаты моделирования при изменении начального значения переменной состояния V на 0,0001 мВ, переходный процесс выделен черным цветом.
Как видно из графиков рис. 6, малое изменение начальных условий привело к значительному изменению длительности переходного процесса в случае модели RK4, в то же время поведение дискретной модели на базе метода DOPRI восьмого порядка не продемонстрировало такой чувствительности. Данный результат указывает на качественную разницу в колебательных режимах этих дискретных моделей. Режим модели RK4 можно охарактеризовать как переходно-хаотический, режим эталонной модели DOPRI – гармонический.
Таким образом, численные эффекты при моделировании нейрона ХХ могут проявляться не только при использовании методов численного интегрирования низких порядков, как было показано в предыдущих экспериментах, но также и при применении распространенного метода Рунге-Кутты четвертого порядка.
Рис. 6. Колебательные режимы RK4 и DOPRI во временной области
Заключение
В данной статье описаны некоторые нежелательные численные эффекты, которые присущи дискретным моделям биологических нейронов. В экспериментах по резонансной генерации потенциалов действия модели нейрона ХХ на основе явного метода Эйлера первого порядка демонстрируют избыточную возбудимость. В ходе экспериментов по поиску хаотических колебаний были обнаружены эффекты, связанные с применением методов интегрирования высоких порядков. Первые результаты по идентификации хаотических режимов модели ХХ при периодическом воздействии были опубликованы в работах [8-10]. Исследования в данной области не потеряли актуальности в наши дни. Так, в работе [11] было изучено поведение классической модели нейрона ХХ при использовании хаотического входного сигнала. В настоящем исследовании особое внимание было уделено явлению переходного хаоса, которое часто проявляется как эффект дискретизации, даже при использовании методов четвертого порядка в случае моделирования нейрона ХХ. Ранее мы исследовали проявления хаотических переходных процессов в модели связанных мемристивными синапсами нейронов ХХ [12], а также в RLC-шунтированной модели джозефсоновского контакта [13].
Результаты настоящей работы могут быть полезны не только при компьютерном моделировании нейронных сетей в инструментальных средах, но также и при цифровой реализации нелинейных динамических систем, описывающих импульсную активность нейронов, к примеру на основе технологии ЦПОС.
1. Jiruska P., De Curtis M., Jefferys J.G., Schevon C.A., Schiff S.J., Schindler K. Synchronization and desynchronization in epilepsy: controversies and hypotheses // The Journal of physiology. 2013. Т. 591, № 4. Р. 787-797.
2. Uhlhaas P.J., Singer W. Abnormal neural oscillations and synchrony in schizophrenia // Nature reviews neuroscience. 2010. Т. 11, № 2. Р. 100.
3. Hodgkin A.L., Huxley A.F. A quantitative description of membrane current and its application to conduction and excitation in nerve // The Journal of physiology. 1952. Т. 117, № 4. Р. 500-544.
4. Бутусов Д.Н., Островский В.Ю., Красильников А.В. Моделирование нелинейных динамических систем параллельными численными методами интегрирования // Фундаментальные исследования. 2014. Т. 9, № 12.
5. Бутусов Д.Н., Каримов А.И., Каримов Т.И., Долгушин Г.К. Семейство аппаратно-ориентированных методов численного интегрирования // Современные проблемы науки и образования. 2014. № 4. С. 138-138.
6. Izhikevich E.M., Moehlis J. Dynamical Systems in Neuroscience: The geometry of excitability and bursting // SIAM review. 2008. Т. 50, № 2. Р. 397.
7. Lai Y.C., Tél T. Transient chaos: complex dynamics on finite time scales. Springer Science & Business Media, 2011. Т. 173.
8. Aihara K., Matsumoto G., Ikegaya Y. Periodic and non-periodic responses of a periodically forced Hodgkin-Huxley oscillator // Journal of theoretical biology. 1984. Т. 109, № 2. Р. 249-269.
9. Guevara M.R., Glass L., Mackey M.C., Shrier A. Chaos in neurobiology // IEEE Transactions on Systems, Man and Cybernetics. 1983. № 5. Р. 790-798.
10. Aihara K. Chaotic oscillations and bifur-cations in squid giant axons // Chaos. 1986. Р. 257-269.
11. Baysal V., Saraç Z., Yilmaz E. Chaotic resonance in Hodgkin-Huxley neuron // Nonlinear Dynamics. 2019. Р. 1-11.
12. Ostrovskii V.Y., Butusov D.N., Belkin D.A., Okoli G. Studying the dynamics of memristive synapses in spiking neuromorphic systems // Conference of Russian Young Researchers in Electrical and Electronic Engineering (EIConRus). IEEE, 2018. Р. 209-214.
13. Ostrovskii V.Y., Karimov A.I., Rybin V.G., Kopets E.E., Butusov D.N. Comparing the Finite-Difference Schemes in the Simulation of Shunted Josephson Junctions // 23rd Conference of Open Innovations Association (FRUCT). IEEE, 2018. Р. 300-305.
