CC BY
18
Вестник КРАУНЦ. Физ.-мат. науки. 2017. № 1(17). C. 33-43. ISSN 2079-6641
DOI: 10.18454/2079-6641-2017-17-1-33-43 МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 517.938
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ НЕРВНОГО ИМПУЛЬСА С УЧЕТОМ ЭРЕДИТАРНОСТИ *
О. Д. Липко
Камчатский государственный университет имени Витуса Беринга, 683032, г. Петропавловск-Камчатский, ул. Пограничная, 4 E-mail: lipko_95@list.ru
В работе предложена математическая модель распространения нервного импульса ФитцХью-Нагумо, которая учитывает эффект эредитарности. Эта эредитарная модель описывается интегро-дифференциальным уравнением со степенным ядром - функцией памяти. Алгоритм численного решения этой модели, реализован в компьютерной программе в среде символьной математики Maple. С помощью этой программы были построены расчетные кривые - осциллограммы, а также фазовые траектории в зависимости от различных значений управляющих параметров.
Ключевые слова: эредитарность, модель ФитцХью-Нагумо, конечно-разностная схема.
© Липко О. Д., 2017
MATHEMATICAL MODELING
MSC 34A08
MATHEMATICAL MODEL OF PROPAGATION OF NERVE IMPULSES WITH
REGARD HEREDITARITY
O.D. Lipko
Vitus Bering Kamchatka State University, 683032, Petropavlovsk-Kamchatsky, Pogranichnaya st., 4, Russia E-mail: lipko_95@list.ru
A mathematical model of the propagation of the nervous pulse of FitzHugh-Nagumo is proposed, which takes into account the effect of heredity. This hereditary model is described by an integro-differential equation with a power kernel - a function of memory. The algorithm for the numerical solution of this model is implemented in a computer program in the environment of symbolic mathematics Maple. With the help of this program, calculated curves - oscillograms, and also phase trajectories were constructed depending on various values of control parameters.
Keywords: hereditarity Model FitzHugh-Nagumo, finite-difference scheme.
© Lipko O.D., 2017
*Работа выполнена по госзаданию, НИР "Применение дробного исчисления в теории колебательных процессов"№AAAA-A17-117031050058-9
Введение
Развитие теории эредитарных динамических систем началось с работы итальянского математика Вито Вольтерра [1], там же он ввел термин эредитарность для описания эффекта последействия, или памяти, и впервые исследовал эредитарный осциллятор. Математическое описание эредитарного осциллятора представляло собой интегро-дифференциальное уравнение с ядром, которое называется функцией памяти. В дальнейшем исследования эредитарных динамических систем были связаны с выбором функции памяти. В силу, того что различные среды могут обладать фрактальными свойствами, то функцию памяти целесообразно выбрать степенной. Тогда интегро-дифференциальные уравнения можно переписать как дифференциальные уравнения дробных порядков, теория которых достаточно хорошо разработана [2]. В литературе такие уравнения называют фрактальными, они описывают процессы с частичной потерей памяти. Фрактальные динамические системы наиболее полно исследовались в монографиях [3, 4].
В работе обобщена динамическая система ФитцХью-Нагумо, которая была предложена Р. ФитцХью [5] и Дж. Нагумо [6] для описания распространения нервного импульса в мембране. Обобщенная математическая модель содержит уравнение с производными дробных порядков в смысле Герасимова-Капуто и решается с помощью конечно-разностной схемы. Основные результаты работы отражены в статье автора [7]. В настоящей работе также, с помощью компьютерных экспериментов, были исследованы вопросы устойчивости и сходимости конечно-разностной схемы, реализующей численное решение предложенной модели.
Постановка задачи и метод решения
Классическая нелинейная динамическая система ФитцХью-Нагумо (ФХН) согласно работам [5], [6] имеет вид:
x3(t)
x (t)= c • (y(t) - x(t)+ z),
. (t)= (x(t) - a + by(t)j (1)
,y (t )= С ,
где a,b, c - константы, удовлетворяющие условиям 1 — 2b/3 < а < 1, 0 < b < 1, b < c2, x(t) - мембранный потенциал, z - интенсивность раздражителя, в первом приближении константа, которая также может иметь вид прямоугольного импульса или дельта-функции, t е [0, T] - время процесса, T > 0 - время моделирования. Динамическая система (1) может быть записана в виде одного уравнения:
x(t) + cx(t) (x2(t) + p)+ qx3 (t) - a - bz = 0, (2)
где p = b/c2 - 1, q = 1 - b, g = b/3. Для уравнения (2) ставятся начальные условия (П, <р - const):
x(0) = n, x (0) = ф. (3)
Задача (2), (3) является задачей Коши, решение которой исследуется в работе
В настоящей работе мы рассмотрим обобщение задачи Коши (2) и (3), введем в него эредитарность с помощью следующего интегро-дифференциального уравнения:
J К\(г - т)х(т)йт - ) + р) I Кг(г - т)х(т+ дх^) + gx3(t) - а - Ъг = 0, (4) о о
где К^ - т) и Kl(t - т) - функции памяти, характеризующие эредитарность.
Замечание. Заметим, что если функции памяти представляют собой дельта-функцию, то тогда в системе отсутствует эредитарность, а если функции памяти представляют собой функции Хэвисайда, то тогда система обладает полной памятью.
Интерес представляет третий вариант: если функции памяти являются степенными функциями, например,
-т) = Тр-«г -т) = Ь-)-),1 <а<2,0<в <1 (5)
где Г(t) - гамма-функция Эйлера, тогда говорят, что система обладает частичной "потерей памяти"[2].
В дальнейшем будем исследовать эредитарные процессы с частичной "потерей памяти". Подставим функции памяти (5) в интегро-дифференциальное уравнение (4). В результате получим:
t 2 t 1 [ Х(т)йт о(х2^)+ р) [ х(т)^т , . 3-. , п ... 1 47 4 4 7 ' 1 к 7 Ьдх(0 + gx3(t) - а - Ъг = 0. (6)
Г(2 - а)У (t - т)а-1 Г(1 - в) У ^ - т)в
Мы получили интегро-дифференциальное уравнение специального вида. Функции памяти (5) в интегро-дифференциальном уравнении (6) могут быть отличны от степенных функций, что приводит к другим интегро-дифференциальным уравнениям. Если обратиться к определению производной дробного порядка по Герасимова-Капуто, то мы приходим к уравнению:
д0"х(т) - с(х2^) + р)дв(х(т) + дх^) + gx3 (t) - а - Ъг = 0, (7)
где дробные дифференциальные операторы равны:
t t дах(т) = 1 [ х(т^т двх(т) = 1 [ х(т^т «т) = г(2 - а)) ^ - т)а-1, ^^ = Г(1 - в)1 (t - т)в , 00
определенные в смысле Герасимова-Капуто с дробными порядками1 < а < 2,0 < в < 1.
Можно отметить, что в предельном случае уравнение (7) переходит в классическое уравнение ФХН (2), поэтому можно считать, что уравнение (2) является частным случаем уравнения (7). Отметим, что уравнение (7) содержит кубическую нелинейность, характерную для осциллятора Дуффинга [8], а также Ван дер Поля [9].
Интегро-дифференциальное уравнение ФХН (7) будем называть дробным, или фрактальным уравнением, а процесс, которые оно описывает, будем называть фрактальным, или эредитарным.
t
t
Задача Коши (7) и (3) в общем виде не имеет точного решения в силу того, что модельное уравнение является нелинейным, поэтому надо использовать численные методы для ее решения. В качестве численного метода возьмем метод конечно-разностных схем, так как его легко можно реализовать в любой компьютерной среде.
Будем рассматривать равномерную сетку. Для этого разобьем временной интервал [0, Т] на N равных частей. В результате получим равномерную сетку tj = ]т, где шаг т = Т/N, j = 0,••• ,N — 1. Значения искомой функции x(tj) = Xj, будем называть ее сеточной функцией. Аппроксимация дробных операторов уравнения (7) осуществляется следующим образом [3, 10]:
k-1
d0tx(x) - wo ^ ■ £ aj ■ (xk-j+1 - 2xk-j + Xk-j-i), aj = (j + l)2 а - j2 а 1(3 - a) j=0
т—в к—1
дЦгХ(т) - в) ■ I ^ ■ j+l — х^—j), ^ = (j + 1)1—в — ^.
Подставим эти аппроксимации в модельное уравнение (7). Приходим к следующей конечно-разностной схеме:
x1 = ф + тп, k = 0, 1
A + Bc(xi + p)
X2 = ~л-72-Ä ((2A + Bc(x2 + p) - q) ■ xi - x3g - Axo + a + bz), k = 1,
21
Хк+1 = А + Вс[х2 + р)((2А + Вс(х2 + р) — д) ■Хк — Хк8 — Ахк—1 + а + Ьг— (8) к—1 к—1 —Вс(х1 + р) ■ I (Ьj(хк+1 — хк)) — А ■ I ^(хк+1 — 2хк + хк—1)), j=l j=l Т—а т—в
А = —--, В = ———, к = 2, ■■■ , п — 1.
Г(3 — а), Г(2 — в), , ,
Замечание. Заметим, что, как правило, нелинейные динамические системы обладают жесткостью при больших значениях управляющих параметров, что приводит к необходимости уменьшить шаг дискретизации в конечно-разностной схеме. В нашем случае, в силу ограниченности параметров а, Ь, с, жесткость отсутствует, поэтому в уменьшении шага нет необходимости.
Результаты моделирования и их обсуждение
Конечно-разностная схема (8) была реализована в компьютерной программе, в среде символьной математики Maple. Рассмотрим применение конечно-разностной схемы (8) численного решения задачи Коши (2) и (3). Значения параметров a,b,c были взяты из работы [5]. Сначала рассмотрим случай, когда изменяются значения дробных параметров а и в, а потом и значения параметра z. Также мы будем исследовать конечно-разностную схему (8) с помощью метода двойного пересчета на ее сходимость и покажем устойчивость по начальным данным и правой части.
Пример 1. Значения управляющих параметров в задаче Коши (2) и (3): t е [0, T ], T = 100, N = 2000, т = 0.05, a = 0.7, b = 0.8, z = -0.4, c = 3, x(0) = 0.2, x (0) = 0.1 Результаты моделирования приведены на рис.1.
Рис. 1. Осциллограммы, полученные по конечно-разностной схеме (8) при значениях параметров а и в: кривая 1 - а = 1.8,в = 0.8, кривая 2 - а = 1.7,в = 0.9, кривая 3 - а = 1.9,в = 0.9, кривая 4 - а = 1.9,в = 0.7
На рис. 1 приведены осциллограммы, полученные по схеме (8) при различных значениях а и в. Осциллограмма под номером 3 по форме похожа на осциллограмму из работы [5]. При уменьшении значений а и в, изменяется форма осциллограмм (смещение периодичности колебаний), однако амплитуда колебаний остается неизменной, что на фазовой плоскости соответствуют предельным циклам (рис.2).
Рис. 2. Фазовые траектории
Исследуем конечно-разностную схему (8) на сходимость с помощью метода двойного пересчета (правила Рунге) при различных значениях параметров а и в .В силу того, аппроксимация уравнения (7) имеет первый порядок, то для вычисления абсолютной ошибки £ мы можем воспользоваться следующей формулой:
£ = тах(|х-Х2г|),I = 1,...,N,
где х{ - численное решение, полученное по схеме (8) на шаге т, х2г -численное решение, полученное по схеме (8) на шаге т/2.
Для оценки расчетной точности можем использовать соотношение:
р = 1п (|£ |) / 1п (т/2) .
Результаты приведены в табл. 1.
Таблица 1
Исследование схемы (8) при различных значениях а ив
а = 1,8 и ß = 0,8
N т Абсолютная ошибка Порядок точности p
10 1/10 0.0456 1.0307
20 1/20 0.0262 0.9868
40 1/40 0.0141 0.9724
80 1/80 0.0073 0.9688
160 1/160 0.0037 0.9691
320 1/320 0.0019 0.9707
а = 1.7 и ß = 0.9
N т Абсолютная ошибка Порядок точности p
10 1/10 0.0591 0.9443
20 1/20 0.0324 0.9301
40 1/40 0.0171 0.9286
80 1/80 0.0088 0.9319
160 1/160 0.0045 0.9364
320 1/320 0.0023 0.9397
а = 1.9 и ß = 0.9
N т Абсолютная ошибка Порядок точности p
10 1/10 0.0436 1.0457
20 1/20 0.0249 1.0003
40 1/40 0.0134 0.9847
80 1/80 0.0069 0.9799
160 1/160 0.0035 0.9790
320 1/320 0.0018 0.9805
а = 1.9 и ß =0.7
N т Абсолютная ошибка Порядок точности p
10 1/10 0.0346 1.1221
20 1/20 0.0202 1.0575
40 1/40 0.0108 1.0313
80 1/80 0.0056 1.0196
160 1/160 0.0028 1.0138
320 1/320 0.0014 1.0119
Из табл.1 можно сделать вывод о том, что при уменьшении шага т, абсолютная ошибка £ уменьшается, а расчетный порядок точности близок к единицы. Такая экспериментальная сходимость не гарантирует сходимости к истинной функции решения задачи Коши (2) и (3). Поэтому необходимо доказывать теорему сходимости.
Пример 2. Рассмотрим другой случай: зафиксируем значения а и в, и будем изменять значения г при значениях параметров: г е [0, Т], Т = 100,N = 2000, т = 0.05, а = 0.7,Ь = 0.8,г = -0.4, с = 3,х(0) = 0.2,х(0) = 0.1, а = 1.8,в = 0.8 и различных значениях г. На (рис. 3) приведены осциллограммы для этого случая.
Рис. 3. Осциллограммы, полученные по конечно-разностной схеме (8): кривая 1 -г = -0.365, кривая 2 - г = -0.4, кривая 3 - г = -0.6, кривая 4 -г = -0.5
В случае г = -0.365 (кривая 1) мы видим, что колебания затухают, а фазовая траектория (рис. 4) имеет вид закручивающейся спирали. При уменьшении значений параметра г , происходит смещение осциллограмм, но с постоянной амплитудой, что обеспечивает выход фазовых траекторий на предельный цикл (рис. 4).
Рис. 4. Фазовые траектории: кривая 1 - г = -0.365, кривая 2 - г = -0.4, ' кривая 3 -г = -0.6, кривая 4 -г = -0.5
Проведем исследования сходимости конечно-разностной схемы (8) в зависимости от значения параметра г. Для этого воспользуемся методикой из предыдущего примера. Результаты приведены а табл. 2.
Таблица 2
Исследование схемы (8) при различных значениях г
z =-0,365
N т Абсолютная ошибка Порядок точности
10 1/10 0.0482 1.0117
20 1/20 0.0278 0.9710
40 1/40 0.0149 0.9588
80 1/80 0.0077 0.9570
160 1/160 0.0039 0.9587
320 1/320 0.0020 0.9615
z =-0,4
N т Абсолютная ошибка Порядок точности
10 1/10 0.0456 1.0307
20 1/20 0.0262 0.9868
40 1/40 0.0141 0.9724
80 1/80 0.0073 0.9688
160 1/160 0.0037 0.9691
320 1/320 0.0019 0.9707
z =-0,5
N т Абсолютная ошибка Порядок точности
10 1/10 0.03775 1.0938
20 1/20 0.02163 1.0391
40 1/40 0.0116 1.0170
80 1/80 0.0060 1.0074
160 1/160 0.0030 1.0031
320 1/320 0.0015 1.0009
z =-0,6
N т Абсолютная ошибка Порядок точности
10 1/10 0.0295 1.1759
20 1/20 0.0168 1.1067
40 1/40 0.0090 1.0742
80 1/80 0.0046 1.0567
160 1/160 0.0023 1.0462
320 1/320 0.0012 1.0393
Аналогично, как и в предыдущем примере, мы видим, что при уменьшении шага т абсолютная ошибка уменьшается £, а расчетный порядок точности стремиться к единице.
Рассмотрим на примере устойчивость конечно-разностной схемы по начальным данным и правой части при следующих значениях управлющих параметров: Т = 1,N = 683,с = 3,а0 = 0.7,г = -0.4,Ь0 = 0.8, а = 1.8,в = 0.8,х(0) = X(0) = 0.2. Для этого сначала добавим малую величину £ = 10-5 в начальное условие х(0), а потом в правую часть уравнения (7). Устойчивость по начальным данным или по правой части будет определяться малым изменением решения задачи Коши (2) и (3) на порядок величины £. В противном случае, решение задачи Коши (2) и (3) будет неустойчивым. Результаты исследования приведены в табл. 3.
Таблица 3
Устойчивость по начальным данным и правой части для схемы (8)
x(0)+ £
т £ 0 £ S
1/500 0.0000100193 0.0000100000 0.0000000193
1/530 0.0000100380 0.0000100000 0.0000000380
1/600 0.0000100543 0.0000100000 0.0000000543
1/685 0.0000100608 0.0000100000 0.0000000608
1/720 0.0000101316 0.0000100000 0.0000001316
f + £
т £ 0 £ S
1/500 0.0000144736 0.0000100000 0.0000044736
1/530 0.0000128115 0.0000100000 0.0000028115
1/600 0.0000122548 0.0000100000 0.0000022548
1/685 0.0000102220 0.0000100000 0.0000002220
1/720 0.0000101831 0.0000100000 0.0000001831
Из табл. 3 мы видим, что для этого примера, имеет место устойчивость по начальным данным и правой части, так как разность 8 между возмущенным и невозмущенным решениями имеет порядок величины е. Конечно, для наиболее полной картины, необходимо доказать теорему об устойчивости конечно-разностной схемы (8). Однако мы в работе показали, что явную-конечно разностную схему можно применять (8) для решения задачи Коши (2) и (3).
Заключение
В работе был предложен и исследован эредитарный нелинейный осциллятор ФитцХью-Нагумо. С помощью теории конечно-разностных схем получено численное решение задачи Коши, построены осциллограммы и фазовые траектории. Показано, что параметры а и в приводят к смещению колебаний осциллятора, но при этом сохраняется постоянство амплитуды, а также изменяется форма фазовых траекторий, которые выходят на предельный цикл. При изменении параметра г колебания могут быть затухающими, а фазовая траектория для соответствующей точки покоя будет являться устойчивым фокусом.
Введение дополнительных управляющих параметров а и в , чтобы более гибко моделировать колебательный режим, дает дополнительную параметризацию сигнала. Дальнейший интерес в исследовании эредитарного осциллятора ФитцХью-Нагумо может заключаться в исследовании на устойчивость точек покоя по аналогии с работой [9], а также дальнейшее обобщение, связанное с введением функций и [11]. Другое направление исследований связано с качественными свойствами конечно-разностной схемы (8) - устойчивостью и сходимостью [10].
Автор выражает благодарность научному руководителю, к.ф.-м.н., Р.И. Паровику за ценные советы и замечания по содержанию данной научной статьи.
Список литературы
[1] Volterra V., "Sur les 'equations int'egro-diff'erentielles et leurs applications", Acta Mathematica, 35:1 (1912), 295-356.
[2] Учайкин В. В., Метод дробных производных, Артишок, Ульяновск, 2008, 512 с. [Uchajkin V. V., Metod drobnyh proizvodnyh, Artishok, Ul'janovsk, 2008, 512 ].
[3] Паровик Р. И., Математическое моделирование линейных эредитарных осцилляторов, КамГУ им. Витуса Беринга, Петропавловск-Камчатский, 2015, 178 с. [Parovik R.
I., Matematicheskoe modelirovanie linejnyh jereditarnyh oscilljatorov, KamGU im. Vitusa Beringa, Petropavlovsk-Kamchatskij, 2015, 178 ].
[4] Petras I., Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation, Springer, Beijing and Springer-Verlag Berlin Heidelberg, 2011, 218 с.
[5] FitzHugh R., "Impulses and physiological states in theoretical models of nerve membrane", Biophysical Journal, 1 (1961), 446-446.
[6] Nagumo J., Arimoto S., Yoshizawa S., "An active pulse transmission line simulating nerve axon", Proc. IRE, 50 (1962), 2061-2070.
[7] Липко О. Д., "Эредитарное модельное уравнение ФитцХью-Нагумо", Международный студенческий научный вестник, 2017, №2, 43-43. [Lipko O. D., "Jereditarnoe model'noe uravnenie FitcH'ju-Nagumo", Mezhdunarodnyj studencheskij nauchnyj vestnik, 2017, №2, 43-43 https://www.eduherald.ru/ru/article/view?id=16890 (дата обращения: 22.04.2017)].
[8] Паровик Р. И., "Математическое моделирование нелокальной колебательной системы Дуффинга с фрактальным трением", Вестник КРАУНЦ. Физико-математические науки, 2015, №1(10), 18-24. [Parovik R.I. Mathematical modeling of nonlocal oscillatory Duffing system with fractal friction. Bulletin KRASEC. Physical and Mathematical Sciences. 2015. vol. 10. no 1. С. 16-21. ].
[9] Паровик Р. И., "Об исследовании устойчивости эредитарного осциллятора Ван дер Поля", Фундаментальные исследования, 2016, №3(2), 283-287. [Parovik R. I., "Ob issledovanii ustojchivosti jereditarnogo oscilljatora Van der Polja", Fundamental'nye issledovanija, 2016, №3(2), 283-287 ].
[10] Parovik R. I., "Explicit finite-difference scheme for the numerical solution of the model equation of nonlinear hereditary oscillator with variable order fractional derivatives", Archives of Control Sciences, 26:3 (2016), 429-435.
[11] Паровик Р. И., "Конечно-разностные схемы для фрактального осциллятора с переменными дробными порядками", Вестник КРАУНЦ. Физико-математические науки, 2015, №2(11), 88-85. [Parovik R.I. Finite-difference schemes for fractal oscillator with a variable fractional order. Bulletin KRASEC. Physical and Mathematical Sciences. 2015. vol.
II. no 2. С. 85-92 ].
Список литературы (ГОСТ)
[1] Volterra V. Sur les 'equations int'egro-diff'erentielles et leurs applications // Acta Mathematica. 1912. vol. 35. issue 1. pp. 295-356.
[2] Учайкин В. В. Метод дробных производных. Ульяновск: Артишок, 2008. 512 c.
[3] Паровик Р. И. Математическое моделирование линейных эредитарных осцилляторов. Петропавловск-Камчатский. КамГУ им. Витуса Беринга. 2015. 178 c.
[4] Petras I. Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation. Beijing and Springer-Verlag Berlin Heidelberg. Springer, 2011. 218 p.
[5] FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane // Biophysical Journal. 1961. vol. 1. pp. 446-446.
[6] Nagumo J., Arimoto S., Yoshizawa S. An active pulse transmission line simulating nerve axon. Proc. IRE. 1962. vol. 50. pp. 2061-2070.
[7] Липко О. Д. Эредитарное модельное уравнение ФитцХью-Нагумо // Международный студенческий научный вестник. 2017. № 2. С. 43-43. url: https://www.eduherald.ru/ru/article/view?id=16890 (дата обращения: 22.04.2017)
[8] Паровик Р. И. Математическое моделирование нелокальной колебательной системы Дуффинга с фрактальным трением // Вестник КРАУНЦ. Физико-математические науки. 2015. №1(10). С.18-24.
[9] Паровик Р. И. Об исследовании устойчивости эредитарного осциллятора Ван дер Поля // Фундаментальные исследования. 2016. № 3(2). C. 283-287.
[10] Parovik R. I. Explicit finite-difference scheme for the numerical solution of the model equation of nonlinear hereditary oscillator with variable order fractional derivatives. Archives of Control Sciences. 2016. vol. 26. issue 3. pp. 429-435.
[11] Паровик Р. И. Конечно-разностные схемы для фрактального осциллятора с переменными дробными порядками // Вестник КРАУНЦ. Физико-математические науки. 2015. №2(11). С. 88-85.
Для цитирования: Липко О. Д. Математическая модель распространения нервного импульса с учетом эредитарности // Вестник КРАУНЦ. Физ.-мат. науки. 2017. № 1(17). C. 33-43. DOI: 10.18454/2079-6641-2017-17-1-33-43
For citation: Lipko O. D. Mathematical model of propagation of nerve impulses with regard hereditarity, Vestnik KRAUNC. Fiz.-mat. nauki. 2017, 17: 1, 33-43. DOI: 10.18454/2079-66412017-17-1-33-43
Поступила в редакцию / Original article submitted: 22.03.2017
