Эффекты дискретизации при численном исследовании модели нейрона Ходжкина - Хаксли

УДК 004.942

DOI: 10.30987/1999-8775-2019-2019-12-94-101

В.Ю. Островский, Д.Н. Бутусов, А.И. Каримов, В.С. Андреев

ЭФФЕКТЫ ДИСКРЕТИЗАЦИИ ПРИ ЧИСЛЕННОМ ИССЛЕДОВАНИИ МОДЕЛИ НЕЙРОНА ХОДЖКИНА - ХАКСЛИ

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

влияние эффектов дискретизации на динамическое поведение биологических моделей нейронов, получены инструменты поиска таких эффектов.

Ключевые слова: модель нейрона, эффекты дискретизации, нелинейная система, численное интегрирование, детерминированный хаос, хаотический переходный процесс.

V.Yu. Ostrovsky, D.N. Butusov, A.I. Karimov, V.S. Andreev

DISCRETIZATION EFFECTS DURING NUMERICAL INVESTIGATION OF HODGKIN-HUXLEY NEURON MODEL

Computer design is a valuable tool in the course of designing neuro-morphic systems. In particular it allows investigating basic mechanisms of neuron pulse activities in networks. For computer modeling it is necessary to digitize a continuous model of the system by means of the application of discrete operators able to keep basic properties of a prototype. But the accuracy of discrete models may decrease because of negative effects caused by the type of the method used, by a discretization pitch and errors in rounding off.

This fact is significant for the analysis of nonlinear systems to which belong the models of biological neurons. As a possible solution of the problem may be the development of specialized tools for the analysis of dynamic systems with the focus upon numerical methods used. In this paper by the example of the modeling of the neuron described by Hodgkin-Huxley classical equations there is considered a set of wide-

Введение

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

spread methods for ODU solution. In the course of investigations there are shown possible negative consequences of incorrect use of some discrete operators. In the paper the results of two sets of computer experiments are presented. The first ones determine the limitations for the practical use of the methods of the first accuracy order during modeling neurons in the mode of resonance generation of action potentials. The second ones show discretization effects connected with chaotic modes of neurons functioning: incorrect behavior of discrete models which is manifested in the emergence of chaotic transition processes. The investigation results may be used at the formation of modeling tool packages both, non-linear dynamic systems in the whole, and neuro-morphic systems in particular.

Key words: neuron model, discretization effects, non-linear system, numerical intergration, deterministic chaos, chaotic transition process.

взгляд на природу импульсной активности нейронов направлен на теорию динамических систем. Ходжкин и Хаксли были первыми, кто изучил динамику нейронов в таком понимании, представив модель мембраны аксона нейрона в 1952 г. [3]. До сих пор модель ионных каналов Ходжкина -Хаксли (ХХ) является наиболее распространенным кинетическим представлением чувствительной к напряжению проводимости в нейронах. Практически все компактные математические модели, демонстрирующие хорошее соответствие экспериментальным данным, основаны на классической модели ХХ.

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

В основе компьютерного моделирования динамических систем лежит применение дискретных операторов, которые предположительно способны сохранить

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

Модель нейрона Ходжкина - Хаксли

Система уравнений Ходжкина - Хаксли [3] представляет собой классическую феноменологическую модель нейрона, устанавливающую связь между внутренними ионными токами и динамическим поведением электрически возбудимой мембраны. Применяя экспериментальный метод фиксации напряжения на гигантском аксоне кальмара, Ходжкин и Хаксли определили, что в клеточной мембране протекают три основных тока: открываемый напряжением устойчивый ток К+, открываемый напряжением устойчивый ток №+ и постоянный ток утечки, переносимый ионами С1-. Эквивалентная электрическая схема для классической модели ХХ показана на рис. 1. Эта схема включает емкость мембраны С, сопротивление утечки Яи две чувствительные к напряжению нелинейные проводимости Сыа и Gк, описывающие ионные каналы натрия и калия, а также три батареи Еыа, Ек и Еь, соединен-

ные последовательно резистивным элементам.

Рис. 1. Эквивалентная электрическая схема модели нейрона ХХ

Система уравнений ХХ для эквивалентной схемы имеет вид

СУ = 1 ~Окп\У-Ек)-0Ыатък{У-ЕЫа)-Оь{У-Еь)-п = ап{У){\-п)-Рп{У)п-т = ат(У)(\-т)-/Зт(У)т-

где

^ (V) = 0.01

10 - V

exp((10 - V) /10) -1

-V

Рп (V) = 0.125 exp | — |;

О, (V) = 0.1

25 -V

exp((25 - V )/10) -1'

Рт (V) = 4 exp

а,(V) = 0.07 exp I |;

20

Ph (V) = 0.1

1

exp((30 - V )/10) +1

В уравнениях системы (1) переменная п моделирует четыре активационных затвора (следовательно, множитель п4) для тока К+, переменная m представляет три

активационных затвора (множитель m3), h является переменной одного инактиваци-онного затвора для тока №+.

Наблюдаемые численные эффекты

Компьютерное моделирование динамической системы (1) выполнено в среде N1 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, Gk = 36 мСм/см2, Gnu = 120 мСм/см2,

Gl = 0,3 мСм/см2, Ek = -12 мВ, Eua = 115 мВ, El = 10,6 мВ. Начальные условия при

моделировании системы (1) - [0,31; 0,05; 0,59; 0].

На рис. 3 представлена бифуркационная диаграмма по параметру Т при т = 5,5 мс. В интервалах 11-12 мс и 16-18 мс дискретная модель, полученная явным ме-

тодом Эйлера с шагом h = 0,05 мс, генерирует потенциалы действия, что является качественным отклонением от прочих моделей на диаграмме, находящихся в режиме подпороговых колебаний.

110-

100" 90-

80"

70-

Ш 60-

¡7 50-

40-

30-

20" ю-

0-

№ И = 0,05

-1 в  ЕМР 7] ООРИ И - 0,005 И = 0,05 И = 0,05

10 11

12 13

14 15

16

17 18 19 20

21

22 23

24 25 26 27 28 29

30

Т [мс]

Рис. 3. Бифуркационная диаграмма для импульсного входного сигнала

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

потенциалов действия для рассматриваемых дискретных операторов. Горизонтальной линией на рис. 4а отмечена строка карты, построенная по бифуркационной диаграмме на рис. 3.

Рис. 4. Динамические карты для импульсного входного сигнала: а - разница между методом первого порядка и эталоном; б - разница между методом второго порядка и эталоном

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

Вторая серия вычислительных экспериментов была проведена для случая, когда динамическая система (1) переходит в хаотический режим генерации потенциалов действия в ответ на сигнал I = ¡0 + А зт(2лю^). Параметры системы (1) остались прежними, за исключением Ема = 120 мВ. Начальные условия установлены на [0,31; 0,05; 0,59; 0,001].

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

сительно их продолжительности I и обычно применяются при изучении статистических особенностей нейронов в ответ на зашумленный входной сигнал. На рис. 5 представлен ряд гистограмм (150 интервалов) для дискретных моделей в ответ на входной ток ¡0 = 6,22 мкА с добавлением 70 Гц синусоидальной компоненты амплитудой А = 0,6 мкА. Время моделирования установлено от 1000 до 3000 мс.

Из рис. 5 видно, что модели на основе EMP, RK4 и DOPRI предоставляют качественно схожие решения, где разброс интервалов свидетельствует о хаотическом режиме на всем временном промежутке моделирования. Дискретная модель, полученная методом EE, на данном промежутке оказалась в периодическом режиме. Анализ данной модели во временной области выявил наличие хаотического переходного процесса, длительность которого чувствительно зависит от начальных условий. Это явление подробно описано в работе [7].

Рис. 5. Межспайковые интервальные гистограммы

Похожий эффект был также найден для модели, полученной методом четвертого порядка ЕК4 при параметрах входного сигнала ¡0 = 6,48 мкА, А = 0,56 мкА и ю = 70 Гц. На рис. 6 показаны результаты моделирования при изменении начального значения переменной состояния V на 0,0001 мВ, переходный процесс выделен черным цветом.

Как видно из графиков рис. 6, малое изменение начальных условий привело к значительному изменению длительности переходного процесса в случае модели

КК4, в то же время поведение дискретной модели на базе метода ПОРЕ1 восьмого порядка не продемонстрировало такой чувствительности. Данный результат указывает на качественную разницу в колебательных режимах этих дискретных моделей. Режим модели ЕК4 можно охарактеризовать как переходно-хаотический, режим эталонной модели ПОРЕ1 - гармонический.

Таким образом, численные эффекты при моделировании нейрона ХХ могут проявляться не только при использовании

методов численного интегрирования низких порядков, как было показано в предыдущих экспериментах, но также и при

применении распространенного метода Рунге-Кутты четвертого порядка.

Рис. 6. Колебательные режимы ЯК4 и БОРЫ во временной области

Заключение

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

сигнала. В настоящем исследовании особое внимание было уделено явлению переходного хаоса, которое часто проявляется как эффект дискретизации, даже при использовании методов четвертого порядка в случае моделирования нейрона ХХ. Ранее мы исследовали проявления хаотических переходных процессов в модели связанных мемристивными синапсами нейронов ХХ [12], а также в Л£С-шунтированной модели джозефсоновского контакта [13].

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

Д.Н. Бутусов, В.Ю. Островский получили поддержку Российского фонда фундаментальных исследований в рамках проекта 19-07-00496/19 «Основы исследовательского проектирования мемристивных систем». А.И. Каримов получил финансирование по гранту Президента для молодых ученых - кандидатов наук МК-811.2019.1.

СПИСОК ЛИТЕРАТУРЫ

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.

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. Butusov D.N., Ostrovsky V.Yu., Krasilnikov A.V. Non-linear dynamic systems modeling by parallel numerical methods of integration // Fundamental Investigations. 2014. Vol.9., No.12.

5. Butusov D.N., Karimov A.I., Karimov T.I., Dol-gushin G.K. Family of hardware oriented methods numerical integration // Current Problems of Science and Education. 2014. No.4. pp. 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. T. 109, № 2. P. 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. P. 790798.

10. Aihara K. Chaotic oscillations and bifur-cations in squid giant axons // Chaos. 1986. P. 257-269.

11. Baysal V., Sarag Z., Yilmaz E. Chaotic resonance in Hodgkin-Huxley neuron // Nonlinear Dynamics. 2019. P. 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. P. 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. P. 300-305.

8. Aihara K., Matsumoto G., Ikegaya Y. Periodic and non-periodic responses of a periodically forced Hodgkin-Huxley oscillator // Journal of theoretical biology. 1984. T. 109, № 2. P. 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. P. 790798.

10. Aihara K. Chaotic oscillations and bifur-cations in squid giant axons // Chaos. 1986. P. 257-269.

11. Baysal V., Sarag Z., Yilmaz E. Chaotic resonance in Hodgkin-Huxley neuron // Nonlinear Dynamics. 2019. P. 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. P. 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. P. 300-305.

Ссылка для цитирования:

Островский В.Ю., Бутусов Д.Н., Каримов А.И., Андреев В.С. Эффекты дискретизации при численном исследовании модели нейрона Ходжкина - Хаксли //Вестник Брянского государственного технического университета. 2019. № 12. С. 94-101. DOI: 10.30987/1999-8775-2019-2019-12-94-101.

Статья поступила в редакцию 21.11.19. Рецензент: к.т.н., доцент Брянского государственного

технического университета, Рытов М.Ю.,

член редсовета журнала «Вестник БГТУ». Статья принята к публикации 29. 11. 19.

Сведения об авторах:

Островский Валерий Юрьевич, ассистент кафедры систем автоматизированного проектирования Санкт-Петербургского государственного электротехнического университета «ЛЭТИ» им. В.И. Ульянова (Ленина), е-шаП: ууоБйоу-skii@etu.ru.

Бутусов Денис Николаевич, к.т.н., доцент кафедры систем автоматизированного проектирования Санкт-Петербургского государственного электротехнического университета «ЛЭТИ» им. В.И. Ульянова (Ленина), е-шаЛ: dnbutusov@etu.ru.

Каримов Артур Искандарович, к.т.н., доцент кафедры систем автоматизированного проектирования Санкт-Петербургского государственного электротехнического университета «ЛЭТИ» им. В.И. Ульянова (Ленина), е-шаЛ: aikarimov@etu.ru. Андреев Валерий Сергеевич, к.т.н., доцент кафедры систем автоматизированного проектирования Санкт-Петербургского государственного электротехнического университета «ЛЭТИ» им. В.И. Ульянова (Ленина), е-шаП: vsandreev@etu.ru.

Ostrovsky Valery Yurievich, Assistant of the Dep. "CAD Systems", Uliyanov(Lenin) State Electro-Technical University "LETI" of Saint-Petersburg, email: vyostrovskii@etu.ru.

Butusov Denis Nikolaevich, Can. Sc. Tech., Assistant Prof. of the Dep. "CAD Systems", Uliyanov(Lenin) State Electro-Technical University "LETI" of Saint-Petersburg, e-mail: dnbutusov@etu.ru.

Karimov Arthur Iskandarovich, Can. Sc. Tech., Assistant Prof. of the Dep. "CAD Systems", Uli-yanov(Lenin) State Electro-Technical University "LETI" of Saint-Petersburg, e-mail: aikarimov@etu.ru. Andreev Valery Sergeevich, Can. Sc. Tech., Assistant Prof. of the Dep. "CAD Systems", Uliyanov(Lenin) State Electro-Technical University "LETI" of Saint-Petersburg, e-mail: vsandreev@etu.ru.