CC BY
7
Нелинейная
^^^^^^^^ динамика и нейронаука
Известия высших учебных заведений. Прикладная нелинейная динамика. 2024. Т. 32, № 3 Izvestiya Vysshikh Uchebnykh Zavedeniy. Applied Nonlinear Dynamics. 2024;32(3)
Научная статья УДК 517.926
DOI: 10.18500/0869-6632-003099 EDN: RJHLMQ
Новый подход к математическому моделированию химических синапсов
Д. С. Глызин, С. Д. ГлызинА.Ю. Колесов
Ярославский государственный университет им. П. Г. Демидова, Россия E-mail: [email protected], [email protected], [email protected] Поступила в редакцию 13.05.2023, принята к публикации 15.01.2024, опубликована онлайн 22.03.2024, опубликована 31.05.2024
Аннотация. Целью данной работы является исследование новой математической модели кольцевой нейронной сети с однонаправленными химическими связями, представляющей собой сингулярно возмущенную систему дифференциально-разностных уравнений с запаздыванием. Методы. С помощью сочетания аналитических и численных методов изучаются вопросы о существовании и устойчивости в этой системе специальных периодических решений — так называемых бегущих волн. Результаты. Предложенные методы позволяют показать, что изучаемая кольцевая система допускает растущее с ростом числа осцилляторов в сети число устойчивых бегущих волн. Заключение. В настоящей статье нами переосмыслен и уточнен предложенный ранее способ математического моделирования химических синапсов. Удалось в полном объеме учесть, с одной стороны, требование вольтерровской структуры соответствующих уравнений и, с другой стороны, гипотезу о насыщающей проводимости. Это позволяет соблюсти принцип единообразия: новая математическая модель строится на тех же принципах, что и предложенная ранее модель электрических синапсов.
Ключевые слова: кольцевая нейронная сеть, химические синапсы, релаксационные циклы, асимптотика, устойчивость.
Благодарности. Исследование выполнено за счет гранта Российского научного фонда № 22-11-00209, https://rscf.ru/ project/22-11-00209/.
Для цитирования: Глызин Д. С., Глызин С. Д., Колесов А. Ю. Новый подход к математическому моделированию химических синапсов//Известия вузов. ПНД. 2024. T. 32, № 3. С. 376-393. DOI: 10.18500/0869-6632-003099. EDN: RJHLMQ
Статья опубликована на условиях Creative Commons Attribution License (CC-BY 4.0).
©Глызин Д. С., Глызин С. Д., Колесов А. Ю., 2024
Article
DOI: 10.18500/0869-6632-003099
A new approach to mathematical modeling of chemical synapses
D.S. Glyzin, S.D. GlyzinM, A. Yu. Kolesov
P. G. Demidov Yaroslavl State University, Russia E-mail: [email protected], El [email protected], [email protected] Received 13.05.2023, accepted 15.01.2024, available online 22.03.2024, published 31.05.2024
Abstract. The purpose of this work is to study a new mathematical model of a ring neural network with unidirectional chemical connections, which is a singularly perturbed system of differential-difference equations with delay. Methods. A combination of analytical and numerical methods is used to study the existence and stability of special periodic solutions in this system, the so-called traveling waves. Results. The proposed methods make it possible to show that the ring system under study allows the number of stable traveling waves to increase with the number of oscillators in the network. Conclusion. In this article, we rethink and refine the previously proposed method of mathematical modeling of chemical synapses. On the one hand, it was possible to fully take into account the requirement of the Volterra structure of the corresponding equations and, on the other hand, the hypothesis of saturating conductivity. This makes it possible to observe the principle of uniformity: the new mathematical model is based on the same principles as the previously proposed model of electrical synapses.
Keywords: circular neural network, chemical synapses, relaxation cycles, asymptotics, stability.
Acknowledgements. This work was supported by the Russian Science Foundation (project No. 22-11-00209), https://rscf.ru/ project/22-11-00209/.
For citation : Glyzin DS, Glyzin SD, Kolesov AYu. A new approach to mathematical modeling of chemical synapses. Izvestiya VUZ. Applied Nonlinear Dynamics. 2024;32(3):376-393. DOI: 10.18500/0869-6632-003099
This is an open access article distributed under the terms of Creative Commons Attribution License (CC-BY 4.0).
1. Описание объекта исследования
Одним из основополагающих принципов построения математических моделей нейронных систем является так называемая гипотеза об эквивалентности. Суть этой гипотезы в том, что мы априори предполагаем эквивалентность биологического нейрона некоторому физическому генератору с сосредоточенными параметрами. В свою очередь, упомянутый генератор моделируется нелинейной системой обыкновенных дифференциальных уравнений или аналогичной системой с запаздыванием. А поскольку колебания мембранного потенциала носят заведомо релаксационный характер, то соответствующая система, как правило, оказывается сингулярно возмущенной.
Упомянутый выше принцип лежит в основе известной модели Ходжкина-Хаксли [1] и многих других математических моделей изолированного нейрона (см. монографию [2] и содержащуюся в ней подробную библиографию). В настоящей работе, придерживаясь гипотезы об эквивалентности, в качестве модели отдельного нейрона мы используем скалярное нелинейное дифференциальное уравнение с запаздыванием вида
и = If (u(t - 1))и (1)
для мембранного потенциала и = u(t) > 0. Здесь параметр X > 0, характеризующий скорость
протекания электрических процессов в нейроне, предполагается большим, точка — дифференцирование по t, а функция f (и) G С2(R+), R+ = {и G R : и ^ 0}, обладает свойствами:
f (0) = 1, f (и) = -а + oQ , uf(u) = o(^) , u2f(u) = O^ (2)
при и ^ +œ>, где а = const > 0. Примером такой функции служит
'(а) = ОТ (3)
Поскольку в предлагаемом ниже способе моделирования химических синапсов уравнение (1) берется за основу, то имеет смысл остановиться вкратце на истории его возникновения и свойствах. В связи с этим обратим внимание на то обстоятельство, что наш подход к моделированию нейронной активности опирается на идеи, принадлежащие Ю. С. Колесову [3] и В. В. Майорову [4], а именно, в монографии [3] описан общий принцип математического моделирования биологических процессов с помощью специальных систем с запаздыванием воль-терровского типа, аналогичных известному уравнению Хатчинсона [5]. Далее, в работе [4] на основе этого принципа и идеи запаздывающей проводимости с насыщением в качестве модели отдельного нейрона было предложено некоторое уравнение с запаздыванием, аналогичное (1). И наконец, в статье [6] после надлежащей доработки упомянутое уравнение приобрело требуемый вид (1), (2).
Следует также отметить, что ранее безотносительно к нейродинамическим приложениям уравнение (1) рассматривалось в работе [7] как одно из возможных обобщений уравнения Хатчинсона. В указанной работе установлено, что при всех X ^ 1 оно допускает экспоненциально орбитально устойчивый цикл u(t, X) > 0, и(0, X) = 1 периода T(X), удовлетворяющий предельным соотношениям:
lim T (X)=To, max |x(î, X) - xo (i)| =o( Л , X ^ +œ, (4)
o 4MT(X) \XJ
где T0 = (1 + a)t0, t0 = 1 + 1/a, x(t, X) = (1/X)ln(«(i, X)), а T0-периодическая функция x0(t) задается равенствами
Ît при 0 < t < 1,
1 -a(t - 1) при 1 < t < to + 1, xo(t + To) = xo(t). (5)
t-To при to + 1 ^t^To,
Наглядное представление о релаксационных свойствах (4) этого цикла дает его график на плоскости (t,и), построенный численно для случая (1), (3) при X = 5, а = 2 (рис. 1), а также график функции (5) (рис. 2).
Обратимся теперь к интересующему нас вопросу о моделировании химических синапсов и напомним, что соответствующая попытка уже предпринималась ранее в статье [8]. Точнее говоря, в [8] был предложен некоторый подход к указанной проблеме, в основе которого лежала подходящим образом модифицированная идея быстрой пороговой модуляции.
Феномен быстрой пороговой модуляции (fast threshold modulation или FTM), впервые описанный в работах [9,10], представляет собой специальный способ связи динамических систем. Характерная особенность этого способа состоит в том, что правые части соответствующих дифференциальных уравнений меняются скачкообразно при переходе некоторых управляющих переменных через свои критические значения. В нейронных системах идея FTM реализуется, как правило, следующим образом.
Рис. 1. График решения x(t) уравнения (1) с функцией (3) при X = 5, а = 2
Fig. 1. Graph of the solution x(t) of the equation (1) with the function f (u) satisfying the formula (3) for X = 5, a = 2
1
\to to - 1 To
0 1 / t
-a
Рис. 2. График функции x0(t) при a = 2 Fig. 2. Graph of the function x0(t) for a = 2
Предположим, что напряжение и = u(t) и сила тока v = v(t) в отдельной нейронной клетке удовлетворяют системе дифференциальных уравнений
eil = f (u,v), v = g(u,v). (6)
Здесь e > 0 — малый параметр, а на правые части f,g <Е С^ наложены стандартные ограничения [11, с. 45-55], обеспечивающие существование устойчивого релаксационного цикла. Типовым примером модели (6) служит известная система ФитцХью-Нагумо [12].
Рассмотрим далее простейшую сеть, состоящую из двух синаптически связанных нейронов. В этом случае согласно выработанным к настоящему времени представлениям (см., например, [13]) соответствующие им электрические переменные (uj ,Vj), j = 1, 2 удовлетворяют системе уравнений
eüi = f (u\,v\) + bs2(u2)(u* - u\), vi = g(ui,vi),
(7)
£il2 = f (U2,V2) + bSi (Ui )(U* - U2), il2 = g(U2 ,Щ).
Здесь b — положительный параметр, характеризующий максимальную проводимость синапса, и* — потенциал покоя (или потенциал Нернста), а функции Sj (uj), j = 1, 2 — постсинаптические проводимости, зависящие от пресинаптических потенциалов Uj.
Следует отметить, что существует несколько различных способов выбора функций Sj (uj), описание которых можно найти в [13]. Мы же, руководствуясь идеей FTM, остановимся на простейшем из них, а именно, будем считать, что
{0 при х < 0,
(8)
1 при х ^ 0,
где и** — порог, начиная с которого одна клетка влияет на другую. Например, если щ < и**, то первый нейрон не действует на второй, а если щ > и**, то действует.
Предположим теперь, что имеет место однонаправленное синаптическое взаимодействие в сети из т, т ^ 2 нейронов, объединенных в кольцо. Если допустить, что каждый отдельно
взятый нейрон моделируется уравнением (1), то мы можем, руководствуясь описанной выше методикой, перейти от (1) к аналогичной (7) системе
й = Xf(uj (t — l))uj + b Sj-\(uj-\)(u* —Uj), j = 1, 2,...,m, (9)
где u0 = um, so = sm, а функции Sj заданы равенствами (8).
Приведенная модель (9) имеет право на существование. Однако, на наш взгляд, в данной ситуации следует отказаться от общепринятых представлений. Первый шаг в этом направлении был сделан в статье [8], где вместо (9) рассматривалась система
Uj = [Xf (uj(t — 1)) + bg(uj-i) ln(u*/u)}uj, j = 1,2,... ,m, uo = um, (10)
в которой b = const > 0, u* = exp(oX), о = const e R, а функция g(u) G C2(R+) такова, что
g(u) > 0 yu> 0, g(0) = 0; g(u) = 1 + О ^ , ug'(u) = О ^ ^
u2g"(u) = О ^^ при u ^
(11)
Мотивы, по которым в [8] мы остановили свой выбор на системе (10), состояли в следующем. Во-первых, при переходе от (9) к (10) общий качественный характер синаптиче-ской связи сохраняется, поскольку в обоих случаях соответствующие связующие слагаемые Ь— из) и Ь\п(и*/и) меняют знак с «+» на «—» при увеличении потенциалов и при прохождении их через критическое значение и*. Во-вторых, в соответствии с методикой работ [3,4] система (10) имеет требуемую вольтерровскую структуру.
Однако со временем стало понятно, что и система (10) нуждается в некоторой доработке. Связано это с тем, что фигурирующая в (10) функция \п(и*/и) не удовлетворяет гипотезе о насыщающей проводимости, согласно которой все входящие в модель нелинейности должны стремиться при и ^ к конечным пределам.
Для того чтобы исправить положение, перейдем от (10) к системе
и з = Х[/(щ (£ — 1)) + Ьд(и^-{)Н(и^/и* )}и,-, 3 = 1, 2,...,т, ио = ит. (12)
Здесь функция Н(и) € С2^+) обладает аналогичными (2), (11) свойствами:
Н(и) > 0 при 0 ^ и < 1, Н(и) < 0 при и> 1, Н(1) = 0, Н'(1) < 1,
Н(0) = 1, Н(и) = —с + О ^^ , иН'(и) = О ^ , и2Н"(и) = О ^^ (13)
при u ^ где с = const > 0.
Свойства (13) гарантируют, что sign h(u) = sign (ln(1/u)) при всех u > 0, а значит, качественный характер связей в моделях (10) и (12) идентичен. Имеет место для h(u) и требуемое свойство насыщения. Однако уместно сделать еще одно нетривиальное допущение, а именно, будем считать, что для каждого уравнения из (12) соответствующий потенциал покоя u*, предполагавшийся постоянным, теперь зависит от времени и совпадает с u,j-i(t). Тем самым приходим к системе
u j = l[f(uj (t — 1)) + bg(uj-i)h(uj/uj-i))uj, j = 1,2,...,m, uo = um, (14)
где, напомним, функции f, g, h обладают свойствами (2), (11), (13) соответственно, b = const > 0, X ^ 1. Данная система и есть интересующая нас новая математическая модель химических синапсов.
Завершая описание объекта исследования, заметим, что предложенная в данной статье модель (14) вполне аналогична рассмотренной в [14] математической модели электрических синапсов, то есть соблюден так называемый принцип единообразия. Именно ради этого принципа и была осуществлена в (12) замена и* на
2. Бегущие волны релейной системы
Характерной особенностью системы (14) является то обстоятельство, что после замен Xj = еlnUj, j = 1, 2,... ,m, e = 1/"k ^ 1 и последующего стремления параметра е к нулю она допускает некоторый предельный объект. Действительно, указанные замены преобразуют ее к виду
±j = F(xj (t — 1),е) + b Q (xj-i,e)H(xj — Xj-l,e), j = 1,2,...,m, (15)
где xo = xm,
F(x,e) = f (exp(Q(x,e) = g (exp (H(x,e) = h (exp (. (16)
Заметим далее, что в силу свойств (2), (11), (13) для функций (16) при любом фиксированном х £ R, х = 0 выполняются предельные равенства
lim F(x, е) = 1 — (а + 1)Н(х), lim Q(х, е) = Н(х),
£—£—(17)
lim Н(х,£) = 1 — (с +1)Н (х), ( )
£—0
где Н(х) — функция из (8).
Соотношения (17) свидетельствуют о том, что при е ^ 0 исходная система (15) переходит в релейную систему
xj = 1 — (а + 1)Н(xj(t — 1)) + ЬН(xj-i)[1 — (с + 1)Н(xj — х—)], 1 < j < т, (18)
где х0 = хт. В свою очередь, наличие предельного объекта (18) существенно облегчает проблему отыскания аттракторов системы (15) и позволяет, в частности, применить к ней общие результаты из [15] о соответствии между грубыми циклами релейной и релаксационной систем. Ниже упомянутые результаты используются для нахождения специальных периодических движений — так называемых бегущих волн.
Согласно общепринятой терминологии бегущей волной с номером к £ N, 1 ^ к ^ т — 1 будем называть специальное периодическое решение системы (15), допускающее представление
Xj = x(t + (j — 1)Д), j = 1,2,...,т. (19)
Здесь Д = const > 0 — некоторый фазовый сдвиг, а функция x(t) — периодическое решение вспомогательного уравнения
х = F (x(t — 1),£) + b д (x(t — Д),е)П(х — x(t — Д),£) (20)
периода Т = тД/к. Аналогичное определение бегущей волны сохраняется и для релейной системы (18), но в этом случае x(t) — периодическое решение периода Т = тД/к вспомогательного уравнения
х = 1 — (а + 1)Н (x(t — 1)) + ЬН (x(t — Д))[1 — (с + 1)Н (х — x(t — Д))]. (21)
Как оказывается, для релейной системы (18) бегущие волны (19) при надлежащем выборе параметров а, Ь, с, к, т могут быть найдены явно. Для того чтобы убедиться в этом, проанализируем уравнение (21), а именно, установим факт существования у него в подходящем диапазоне параметров а, Ь, с, А некоторого периодического решения х = х*(Ъ, А) периода Т* = Т*(Д).
Предполагая выполненным условие А > 1, в фазовом пространстве С [—А, 0} непрерывных при —А ^ £ ^ 0 скалярных функций ф(£) введем в рассмотрение специальное семейство Фт1,т2 (¿) начальных условий для уравнения (21). Это семейство зависит от двух вспомогательных параметров Т1, Т2, удовлетворяющих неравенствам
0 < Т1 <
а + 6+1
Т2,
Т2 < 1,
6 + 1
+ Т2 + 1 < А < Т2 + 1 + а.
(22)
Сами же функции фТ1,Т2 (t) задаются формулой
t + А — т1 + а(т2 — т1) при —А —А + т1,
ФГ1,Т2 № = <
—a(i + А — Т2)
— а + + А — Т2 — 1
(6+1)i
при —А + т1 ^ t ^ —А + т2 + 1, при —А + т2 + 1 ^ t ^ —А + Тз, при —А + т3 ^ t ^ 0,
где
тз = А + -(А — Т2 — 1 — а). Обратим внимание, что в силу соотношений (22)-(24) справедливы неравенства
—А < —А + Т1 < —А + Т2 < —А + Т2 + 1 < —А + тз < 0, фт1 ,т2 (¿) > 0 при — А < —А + Т2, Фт1,т2 (^ ) < 0 при — А + Т2 <г< 0,
(23)
(24)
(25)
а значит, график любой функции фТ1,Т2 (£) из нашего семейства имеет вид, показанный на рис. 3.
Обозначим через х
Т1 ,Т2
(t), t ^ 0 решение уравнения (21) с начальным условием
фт1,т2 (¿), —А ^ £ ^ 0, а через £ = Т > 0 — второй положительный корень уравнения хт1,т2 (¿) = 0 (если он существует). Как будет показано ниже, при некоторых дополнительных ограничениях
на а, , , А и на параметры т1 , т2 функция хт1,т2 (£ + Т), —А ^ £ ^ 0 совпадает с Фтьт2 (¿), где новые параметры Т1, Т2 линейно зависят от Т1, т2. Более того, соответствующее отображение (т1, т2) ^ (т1, т2) имеет единственную неподвижную точку (т1 (А),т2(А)). Что же касается искомого периодического решения х* (£, А) уравнения (21), то само оно и его период Т* (А) задаются соотношениями
x*(t, А) — ХТ1,Т2 (^)|Т1=Т*(А),Т2=Т*(А), Т*(А) — Т |Т1=Т*(А),Т2=Т2(А).
(26)
Рис. 3. График функции фТ1,Т2 (t) из семейства (23) при условиях (22), (24), (25)
Fig. 3. Graph of function фТ1,Т2 (t) from the family (23) under conditions (22), (24), (25)
Для того чтобы реализовать описанную выше программу действий, проинтегрируем уравнение (21) методом шагов, а именно, будем последовательно рассматривать различные промежутки времени и получать для интересующего нас решения хт1,т2 (£) на этих промежутках некоторые явные формулы.
а
а
На первом шаге обратимся к промежутку
а
0 < и, Ь* = Т2. (27)
а +Ь+1
В этом случае в силу неравенств (22), (25) имеем
ь- 1 е [-1, и - 1) с (-Д + Т2,0), х(г - 1) = фтьт2^ - 1) < 0,
а значит, Н(х(Ь - 1)) = 0. Далее, учитывая вытекающее из (22), (27) неравенство %\ < ¿*, приходим к формуле
т1 + а(т2 - Т1) при 0 XI,
х(1 - Д) = фтх,Т2 ^ - Д) = и 1 ( 2 1) " " 1 (28)
I- а^ - т2) при т1 ¿*.
А отсюда и из (25) автоматически следует, что х(Ь - Д) > 0, Н(х(Ъ - Д)) = 1. И, наконец, в случае функции х(Ь) - х(Ъ - Д) мы априори считаем выполненным неравенство
х(Ь) -х(г - Д) < 0 УЬе [0, и). (29)
Принимая во внимание перечисленные факты, приходим к выводу, что на первом шаге рассмотрению подлежит задача Коши
х = Ъ + 1, х\1=0 = фт1;Т2 (0) = 0.
Тем самым на промежутке (27) решение х = хТ1Т2 (Ь) задается равенством
х = (ь + 1)г. (30)
Следует, однако, напомнить, что это равенство выведено в предположении (29). Но из формул (28), (30) нетрудно вывести, что упомянутое априорное условие действительно выполняется. Второй шаг связан с рассмотрением интервала
Ь* < Ь < Ь**, Ь** = Т2. (31)
Отметим сразу, что в точке £ = Ь* решение хТ1Т2 (Ь) доопределяется по непрерывности, то есть хТ1,Т2 (к*) = (Ь + 1)к*. Учитывая это обстоятельство, приходим к выводу, что при £ = Ь* соответствующее неравенство из (29) обращается в строгое равенство. Поэтому уместно предположить, что в точке = * происходит переключение и
х(ь) -х(г - Д) > 0 Уье (г*, и*). (32)
Проанализируем далее знаки функций х(Ъ - 1), х(Ь - Д). Как и в предыдущем случае, для первой из них имеем
ь-1 е (г* - 1, и* -1) с (-Д + Т2,0), х(г - 1) = фтьт2 ^ -1) < 0, н(х(г -1)) = 0. Вторая же на промежутке (31) задается формулой
х(ъ - д) = фтьт2^ - Д) = -а(г - г**), (33)
из которой, в свою очередь, следует, что х^ - Д) > 0, Н(х(Ъ - Д)) = 1.
Приведенные факты свидетельствуют о том, что на интервале (31) имеем дело с задачей
Коши
х = 1 -Ьс, х\г=и =(Ь+1)г *,
а значит, в этом случае для решения х = хТ1 ,Т2 (Ь) справедливо равенство
х = (Ь + 1)и - (Ьс- 1)(г - и). (34)
Объединяя его с (33), приходим к выводу, что
х(г) -х(г - Д) = (а + 1 -Ъс)(г - и) Vг е [г*, и*).
Тем самым априорное требование (32), при котором была получена формула (34), заведомо справедливо при условии (которое всюду ниже считаем выполненным)
Ьс<а + 1. (35)
Добавим еще, что, как обычно, в точке £ = Ь** решение хТ1 ,Т2 (к) доопределяется по непрерывности, а условие (35) гарантирует, что
I л_(&+1)(а + 1 -Ьс)
хТ1,Т2**) — а +1) +1 > (36)
Прежде чем перейти к последующему построению решения хТ1 ,Т2 (к), сделаем одно полезное наблюдение. Поскольку в силу (25) справедливо неравенство х( - Д) < 0 при е ( **, Д), то автоматически Н(х(Ъ - Д)) = 0. Тем самым на промежутке (Ь**, Д) вместо (21) рассмотрению подлежит более простое уравнение
х = 1 - (а + 1)Н(х(г - 1)). (37)
На третьем шаге обратимся к интервалу
и* <г< 1 (38)
и заметим, что согласно (25) при указанных £ справедливы соотношения х(Ь-1)<0, Н(х(Ъ-1)) = 0. А отсюда и из (36), (37) заключаем, что на промежутке времени (38) решение хТ1,Т2 (Ь) задается формулой
х = к * + ■ : ' " t**, (39) а + + 1
а в точке Ь = 1 доопределяется по непрерывности.
Перед рассмотрением очередного промежутка времени положим
1 ( (Ь + 1)(а + 1 -Ьс) \
и** = 1 + -[1 -г* * + ------ и * > 1 (40)
а а + + 1
и будем считать выполненным очередное условие
и** + 1 < Д. (41)
Четвертый шаг связан с интервалом
1 <г< г*** +1, (42)
на котором в силу (41) по-прежнему имеем дело с уравнением (37). Более того, предполагаем априори, что при указанных значениях времени
х(* — 1) > 0. (43)
Принимая во внимание условие (43), из (37), (39), (40) заключаем, что здесь для решения хт1,т2 (£) имеет место равенство
х = — а( — ***). (44)
Далее, объединяя формулу (44) с уже полученными ранее формулами для хт1,т2 (£) на первых трех шагах (см. (30), (34), (39)), приходим к выводу о справедливости на промежутке времени (42) требуемого условия (43). Тем самым равенство (44) обретает законную силу, а в точке £ = Ь*** + 1 по непрерывности получаем
хТ1,Т2 *** + 1) = а. (45)
На пятом шаге обратимся к интервалу
и** + 1 <г < А, (46)
считая, что при указанных справедливо неравенство
х( — 1) < 0. (47)
Далее, учитывая в (37) вытекающее из (47) тождество Н(х(Ъ — 1)) = 0 и дополняя получившееся уравнение начальным условием (45), для хт1,т2 (£) приходим к очередной формуле
х = Ь — Ь*** — а — 1. (48)
Предположим, что выполнено условие
А < Ь*** +а + 1. (49)
Тогда в силу формул (44), (48) автоматически выполняется и априорное требование (47). Тем самым на интервале (46) равенство (48) действительно имеет место и по непрерывности
хТ1 ,Т2 (А) = А — I*** — а — 1 < 0. (50)
Прежде чем перейти к заключительному этапу построения решения хт1,т2 (£), введем в рассмотрение величину
Т = А + *** + а + 1 — А), (51)
предполагая дополнительно, что
А + и* <Т< А + 1. (52)
Подчеркнем, что согласно (51) условие (49) означает, что Т > А. Таким образом, требования (52) более сильные, чем ограничение (49).
На шестом шаге рассмотрим промежуток
А <г ^Т (53)
в предположении о выполнении неравенств
х(£) < 0, х(Ъ — 1) < 0, х(Ъ — А) > 0. (54)
Далее, учитывая информацию (54) в правой части уравнения (21) и опираясь на формулы (50), (51), для отыскания хТ1 ,т2 (£) приходим к задаче Коши
х = Ь + 1, х|4=А = (6 + 1)(А — Т).
Таким образом, в случае (53) решение хт1,т2 (¿) имеет вид
ж = (b + 1)(t — Т).
(55)
Что же касается априорных условий (54), то их справедливость гарантируют неравенства (52), формула (55) и уже полученные ранее явные выражения для хт1,т2 (£) на предыдущих пяти шагах (см. (30), (34), (39), (44), (48)).
Наглядное представление о процессе построения решения хт1,т2 (£) уравнения (21) дает рис. 4, где показан график хТ1 ,т2 (¿) на отрезке 0 ^ £ ^ Т .Из этого графика и из условий (52), в частности, следует, что функция хТ1 ,т2 (£ + Т) имеет на отрезке —А ^ £ ^ 0 ровно три излома (как и исходное начальное условие фт1,т2 (£)). Более того, опираясь на соотношения (39), (44), (48), (55), нетрудно заметить, что хт1,т2(£ + Т) = фт1,т2(¿), —А ^ ^ ^ 0, где
Рис. 4. График хХ1 ,Т2 (t) на отрезке 0 < t < Т Fig. 4. Graph of function жТ1 ,T2 (t) for t G [0, T]
Ti = 1 + A — Т, T2 = A + U** - Т.
(56)
Покажем теперь, что отображение (т1, Т2) ^ (Т1, Т2), порожденное соотношениями (56), имеет единственную неподвижную точку (т*(А), т2(А)). С этой целью подставим в (56) явные формулы для Ь*** и Т (см. (40), (51)). В результате после некоторых преобразований приходим к выводу, что компоненты интересующей нас неподвижной точки задаются равенствами
T* (A) = iTÏ (A + 'a— 1 — \
+
b( a — с(Ь + 1)) a(b + 1)(a + 6 + 1)
5(A),
(57)
т2 (A) =
а + 6 + 1
а(а + 1)(6 + 1) + b(a + cb(b + 1))
(aA — (a + 1)(a — b)).
(58)
Подведем некоторый итог. Обозначим через совокупность всевозможных наборов параметров (а,Ь,с, А), удовлетворяющих условиям (22), (35), (41), (52), в которых т1, т2 заданы равенствами (57), (58). Из проделанных выше построений вытекает следующая
Теорема 1. При любом наборе (а,Ь,с, А) € уравнение (21) имеет периодическое решение х = х*(£,А) периода Т*(А), задающееся формулами (26), (57), (58). Более того, для Т*(А) справедливы соотношения
Т* (A) = GiA + 62, е
i=
6 + 1
1+
a — c(b + 1)
a(a + 1)(6+ 1) + b(a + cb(b + 1))/
G2 =
a + 1 a(b + 1)
a + 1 —
( a — b)b(a — c(b + 1)) a(a + 1)(6 + 1) + b(a + cb(b + 1))
(59)
T
Установленная теорема позволяет получить в явном виде некоторые бегущие волны релейной системы (18). Для того чтобы сделать это, фиксируем произвольно натуральное к : 1 < к < т - 1 и, опираясь на формулы (59), из уравнения Т*(Д) = тД/к найдем
Д*'"я = т/кЪн ■ (60)
Далее, через Q2 обозначим множество наборов ( а, Ь, с, т, к), для которых
т/к - 61 = 0, Д*п'к) > 0, (а,Ь,с, Д)\Д=Д(тЛ) е (61)
где Д— величина (60). Нетрудно заметить, что при условиях (61) система (18) допускает бегущую волну вида (19) с номером к. Эта волна задается равенствами
х, = х{Г'к\1 + и - 1)Д(Г'к)), 3 = 1,2,... ,т, (62)
где х*Г'к){1) = х*(1, Д)|Д=Д(т,*о, и имеет период Т*т'к) = тД*п'к^)/к.
3. Бегущие волны релаксационной системы
Перейдем теперь к вопросу о существовании аналогичной (62) бегущей волны у исходной релаксационной системы (15). Для этого сначала обратимся к соответствующему вспомогательному уравнению (20). Аналогом теоремы 1 здесь является следующее утверждение.
Теорема 2. Предположим, что фигурирующие в (16) функции ¡(и), д(и), к(и) удовлетворяют условиям монотонности
¡'(и) < 0, д'(и) > 0, к'(и) < 0 Vи ^ 0. (63)
Тогда по любому набору параметров (а, Ь,с, Д) е можно указать такое достаточно малое ео > 0, что при указанных а, Ь, с, Д и при всех 0 < е < ео уравнение (20) допускает периодическое решение х*(1 ,е, Д), х*(0,е, Д) = 0 периода Т*(е, Д) с асимптотикой
Т* (е, Д) =Т*(Д)+0(е), тах \х*(Ь,е, Д) - х*(Ъ, Д)\ =0(е), е^ 0. (64)
ге{о,т*(е,Д)]
Здесь х*(1, Д) — периодическое решение уравнения (21) периода Т*(Д), доставляемое теоремой 1.
Поскольку доказательство сформулированной теоремы базируется на общих результатах из [15], то ограничимся лишь описанием общей его схемы, состоящей из трех этапов.
На первом этапе зафиксируем достаточно малые оо, о1 > 0 и введем в рассмотрение множество 5 начальных функций ф(£) е С [-Д - оо, - оо], удовлетворяющих условиям
\ф^) - х* (г, Д)\ < 01 при - Д - Оо -оо, ф(—оо) = -(Ь+1)оо.
Далее, обозначим через хф^,е, Д), £ ^ -оо — решение уравнения (20) с произвольной начальной функцией ф е 5. Опираясь на условия монотонности (63) и аппарат дифференциальных неравенств, удается показать, что второй положительный корень £ = Тф(е, Д) уравнения
Жф(t — Оо,е, А) = —(6+1)оо существует и является простым. Тем самым на S корректно определен оператор последования Пуанкаре
Пе(ф)=Хф(г + Тф(е,А),е,А), —А — оо —оо (65)
со значениями в С [—А — оо, — оо]. Полное доказательство факта существования и простоты второго положительного корня уравнения для определения оператора Пуанкаре П£(ф) приведено, например, в статьях [14] и [16]. Это доказательство с незначительными изменениями технического характера переносится и на рассматриваемый случай. Учитывая, что соответствующие выкладки довольно громоздки, они здесь не приводятся.
На втором этапе введем в рассмотрение решение x(t,е, А) уравнения (20) с начальной функцией x*(t, А), —А — оо ^ t ^ — оо. Асимптотическое интегрирование соответствующей задачи Коши приводит к выводу, что при е ^ 0 имеют место асимптотические равенства
Тф(£, А)|ф=л*(4>А) = Т* (А) + О(е), max \x(t,£, А) — x*(t, А)| =0(е). (66)
Соотношения (66) означают, что функция x = x(t,е, А) является приближенным (с точностью до О( )) периодическим решением уравнения (20).
На третьем этапе рассмотрим линейный оператор
А(е, А) = 9фП£(ф)|ф=ж*(4;А) : Со [—А — оо, —оо] ^ Со [—А — оо, —оо], (67)
где 9фПе(ф) — производная Фреше оператора (65),
Со[—А — оо, —оо] = {ф(£) е С [—А — оо, —оо ] : ф(—оо) = 0}.
Асимптотический анализ отвечающего уравнению (20) линейного уравнения в вариациях на решении x = x(t, е, А) позволяет заключить, что оператор (67) допускает простое вещественное собственное значение X = Хо (е):
1- w ч &2( а — с( Ь + 1))
lim Хо(е ) = ^—г-т-^—т ч < 1. (68)
в^о ок J а(Ь + 1)(а + 6 + 1) v ;
Остальной же спектр этого оператора лежит в шаре {X е C : |Х| ^ т(е)}, где r(e) ^ 0 при е ^ 0.
Подводя итог, отметим, что соотношения (66)-(68) гарантируют применимость к уравнению П£(ф) — ф = 0 в точке (е, ф) = (0, x*(t, А)) теоремы о неявном отображении. Тем самым оператор (65) допускает неподвижную точку ф*(£,е, А) е S, асимптотически близкую к x*(t, А). Ясно также, что в исходном уравнении (20) этой неподвижной точке соответствует периодическое решение, которое после надлежащего сдвига по времени перейдет в искомое периодическое решение x*(t,е, А) с асимптотикой (64).
Из первого равенства (64) очевидным образом следует, что в случае
(а, Ь, с, т, к) е Ü2 (69)
при всех достаточно малых е > 0 уравнение Т*(е, А) = тА/к допускает решение
А*т'к)(е)= А*т'к) +0(е), е^ 0. (70)
В свою очередь, подставляя соотношение (70) в ж* (t,е, А), получаем функцию , е), пе-
риодическую с периодом Т(т'к\е) = тА*п'к^ (е)/к. А это значит, что нами построена бегущая волна системы (15) вида
x, = x{r,k)(t + (j — 1)А{Г'к)(е), S), j = 1, 2,...,т. (71)
В принципе, вопрос об устойчивости цикла (71) может быть решен теоретически (см., например, [8, 17]). Однако получающиеся на этом пути формулы не удается проанализировать аналитическими методами. В связи с этим нами предпринята серия численных экспериментов, состоящая из двух этапов.
На первом этапе фиксируются параметры а = 2, b = 1, с = 2 и при изменении числа осцилляторов т от 10 до 50 выясняется, при каких значениях к для набора параметров (а, Ь, с, т, к) справедливо включение (69). Выбор величин а,Ь,с в значительной степени случаен, единственное требование состоит в том, чтобы при различных значениях количества осцилляторов т, число устойчивых бегущих волн было по возможности большим.
На рис. 5 изображена зависимость величины N(т), определяемой числом подходящих значений к, от параметра 10 ^ т ^ 50. Оказалось, что N(т) растет с ростом т по близкому к линейному закону.
На втором этапе для наборов параметров, которым соответствуют бегущие волны, численно интегрируется релейная система уравнений с запаздыванием (18). Нами используется метод VODE, который проверяет задачу на жесткость и для нежестких задач представляет собой неявный метод Адамса, а для жестких — метод из класса BDF (backward differentiation formulas) (см., [18]). В качестве начальных условий для этой системы выбираются функции (62) при —1 ^ t ^ 0 с добавлением малого случайного возмущения. Напомним, что эти функции известны нам в явном виде. На основании результатов этого интегрирования делается вывод об устойчивости или неустойчивости бегущей волны с номером к.
Вычисления показали, что для числа осцилляторов т от 5 до 20 оказывается устойчивой лишь одна из имеющихся бегущих волн. При т ^ 21 удается получить две устойчивые волны. Большего числа сосуществующих устойчивых волн в рамках данного эксперимента получить
N
m
20
30
40
50
Рис. 5. Число N(т) устойчивых бегущих волн (71) (подходящих значений к) при 10 < гп < 50 Fig. 5. Number N(т) of stable traveling waves (71) (suitable values of k) at 10 < rn < 50
0.5 X у.
0
-0.5
-1.0
-1.5
-2.0 ^ t
19 20 21 22 23 24 25 26 Рис. 6. Зависимость от t первой компоненты xi(t) устойчивой бегущей волны системы (18) при гп = 19 и к = 16 Fig. 6. Dependence on t of the first component xi(t) of a stable traveling wave of the system (18) for rn = 19 and к = 16
не удалось. В то же время численный счет показал, что релейная система (18), а значит, и релаксационная система (15), имеют и другие устойчивые периодические режимы.
На рис. 6 изображена зависимость от £ первой компоненты х\(1) устойчивой бегущей волны системы (18) при т = 19 и к = 16, все остальные компоненты решения этой системы представляют собой ту же самую функцию х\(1) с соответствующим номеру волны сдвигом (см. формулу (71)).
Заключение
В настоящей статье мы переосмыслили и уточнили предложенный в [8] способ математического моделирования химических синапсов. В отличие от модели (10), в новой ее версии (14) в полном объеме учтены требования вольтерровской структуры соответствующих уравнений и гипотеза о насыщающей проводимости. Кроме того, как уже отмечалось выше, соблюден и так называемый принцип единообразия: новая математическая модель (14) аналогична предложенной в [14] модели электрических синапсов.
Завершая рассмотрение модели (14), отметим две нерешенные задачи. Первая из них связана с пороговым значением u = u** из (8). В случае системы (14) мы изначально считали, что аналог этого значения равен единице. Если же это не так, то вместо (14) приходим к системе
Uj = A[f (Uj (t — 1)) + bg(uj-i/u**)h(uj/uj-i)]uj, j = 1,2,... ,m, (72)
где u0 = um, u** = exp(A d), d = const £ R.
Как и в случае (14), системе (72) отвечает аналогичная (18) релейная система
x-j = 1 — (а + 1)Н(xj(t — 1)) + ЬН(xj—i — d)[1 — (с + 1)Н(xj — х—)], 1 < j < m, (73)
где хо = хт. Представляет интерес проблема отыскания устойчивых бегущих волн сначала у предельного объекта (73), а затем и у исходной системы (72). Ясно, что по аналогии с проделанным в данной работе анализом, сформулированная проблема вполне может быть решена посредством надлежащего сочетания аналитических и численных методов.
Другая нерешенная задача состоит в исследовании аттракторов аналогичной (14) полносвязной сети взаимодействующих нейронов. Соответствующая математическая модель имеет вид
iij = А,
f(ui (t — ^ + bY^9(us)h(uj /us)
S=1
S=j
j = 1,2,...,m, (74)
где функции f, д, h те же самые, что и в (14), b = const > 0, X ^ 1. Согласно общим результатам из [16,19] в системе (74) могут существовать как периодические режимы двухкластерной синхронизации, так и бегущие волны. Проблемы отыскания указанных периодических движений пока не решены.
Список литературы
1. Hodgkin A. L., Huxley A. F. A quantitative description of membrane current and its application to conduction and excitation in nerve // Journal of Physiology. 1952. Vol. 117. P. 500-544. DOI: 10.1113/jphysiol.1952.sp004764.
2. Ижикевич Е. МДинамические системы в нейронауке. Геометрия возбудимости и пачечной активности. М.; Ижевск: Ижевский ин-т компьютерных исследований, 2018. 520 с.
3. Колесов А. Ю., Колесов Ю. С. Релаксационные колебания в математических моделях экологии. Тр. МИАН, 199. М.: Наука, 1993. 126 с.
4. Майоров В. В., Мышкин И. Ю. Математическое моделирование нейронной сети на основе уравнений с запаздыванием // Матем. моделирование. 1990. Vol. 2, no. 11. С. 64-76.
5. Hutchinson G. E. Circular causal systems in ecology // Ann. N. Y. Acad. of Sci. 1948. Vol. 50, no. 4. С. 221-246. DOI: 10.1111/j.1749-6632.1948.tb39854.x.
6. Глызин С. Д., Колесов А. Ю., Розов Н. Х.Релаксационные автоколебания в сетях импульсных нейронов // УМН. 2015. Т. 70, №3(423). С. 3-76.
7. Колесов А. Ю., Мищенко Е. Ф., Розов Н. Х. Об одной модификации уравнения Хатчинсона // Ж. вычисл. матем. и матем. физ. 2010. Vol. 50, no. 12. С. 2099-2112.
8. Глызин С. Д., Колесов А. Ю., Розов Н.Х. Об одном способе математического моделирования химических синапсов // Дифференц. уравнения. 2013. Vol. 49, no. 10. С. 1227-1244. DOI: 10.1134/S0374064113100014.
9. Somers D., Kopell N. Rapid synchronization through fast threshold modulation // Biol. Cybern. 1993. Vol. 68. P. 393-407. DOI: 10.1007/BF00198772.
10. Kopell N., Somers D. Anti-phase solutions in relaxation oscillators coupled through excitatory interactions // J. Math. Biol. 1995. Vol. 33. P. 261-280. DOI: 10.1007/BF00169564.
11. Мищенко Е. Ф., Розов Н.Х. Дифференциальные уравнения с малым параметром и релаксационные колебания. М.: Наука, 1975. 248 с.
12. FitzHugh R. A. Impulses and physiological states in theoretical models of nerve membrane // Biophysical J. 1961. Vol. 1. P. 445-466. DOI: 10.1016/S0006-3495(61)86902-6.
13. Terman D., Borisyuk A, Friedman A, Ermentrout B. An Introduction to Dynamical Systems and Neuronal Dynamics // Tutorials in Mathematical Biosciences I. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer-Verlag, 2005. P. 21-68. DOI: 10.1007/978-3-540-31544-5_2.
14. Глызин С. Д., Колесов А. Ю. Об одном способе математического моделирования электрических синапсов // Дифференц. уравнения. 2022. Vol. 58, no. 7. С. 867-881. DOI: 10.31857/S037 4064122070019.
15. Колесов А.Ю., Мищенко Е. Ф., Розов Н.Х. Реле с запаздыванием и его С1 -аппроксимация // В сб.: Динамические системы и смежные вопросы: К 60-летию со дня рождения академика Дмитрия Викторовича Аносова. Тр. МИАН. 1997. Vol. 216. С. 126-153.
16. Глызин С. Д., Колесов А. Ю. Бегущие волны в полносвязных сетях нелинейных осцилляторов //Ж. вычисл. матем. и матем. физ. 2022. Vol. 62, no. 1. С. 71-89. DOI: 10.31857/S00444669 22010070.
17. Глызин С. Д., Колесов А. Ю., Розов Н. Х. Явление буферности в кольцевых генных сетях // ТМФ. 2016. Vol. 187, no. 3. С. 560-579. DOI: 10.4213/tmf9052.
18. Brown P. N, Byrne G. D., and Hindmarsh A. C. VODE: A Variable Coefficient ODE Solver // SIAM J. Sci. Stat. Comput. 1989. Vol. 10, no. 5. P. 1038-1051. DOI: 10.1137/0910062.
19. Глызин С. Д., Колесов А.Ю. Периодические режимы двухкластерной синхронизации в полносвязных сетях нелинейных осцилляторов // ТМФ. 2022. Т. 212, № 2. С. 213-233. DOI: 10.4213/tmf10191.
References
1. Hodgkin AL, Huxley AF. A quantitative description of membrane current and its application to conduction and excitation in nerve. Journal of Physiology. 1952;117:500-544. DOI: 10.1113/ jphysiol.1952.sp004764.
2. Izhikevich Eugene M. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. Cambridge, Massachusetts, USA: The MIT Press; 2007. 464 p.
3. Kolesov AY, Kolesov YS. Relaxational oscillations in mathematical models of ecology. Proc. Steklov Inst. Math. 1995;199:1-126.
4. Maiorov VV, Myshkin IY. Mathematical modeling of a neuron net on the basis of the equation with delays. Math. Models Comput. Simul. 1990;2(11):64-76 (in Russian).
5. Hutchinson GE. Circular causal systems in ecology // Ann. N. Y. Acad. of Sci. 1948;50(4): 221-246. DOI: 10.1111/j.1749-6632.1948.tb39854.x.
6. Glyzin SD, Kolesov AYu, Rozov NKh. Self-excited relaxation oscillations in networks of impulse neurons. Russian Math. Surveys. 2015;70(3):383-452. DOI: 10.1070/RM2015v070n03ABEH 004951.
7. Kolesov AYu, Mishchenko EF, Rozov NKh. A modification of Hutchinson's equation. Comput. Math. Math. Phys. 2010;50(12):1990-2002. DOI: 10.1134/S0965542510120031.
8. Glyzin SD, Kolesov AYu, Rozov NKh. On a method for mathematical modeling of chemical synapses. Differential Equations. 2015;70(3):383-452. DOI: 10.1134/S0012266113100017.
9. Somers D, Kopell N. Rapid synchronization through fast threshold modulation. Biol. Cybern. 1993;68:393-407. DOI: 10.1007/BF00198772.
10. Kopell N, Somers D. Anti-phase solutions in relaxation oscillators coupled through excitatory interactions. J. Math. Biol. 1995;33:261-280. DOI: 10.1007/BF00169564.
11. Mishchenko EF, Rozov NKh. Differential Equations with Small Parameters and Relaxation Oscillations. New York, NY: Springer; 1980. 228 p. DOI: 10.1007/978-1-4615-9047-7.
12. FitzHugh RA. Impulses and physiological states in theoretical models of nerve membrane. Biophysical J. 1961;1(6):445-466. DOI: 10.1016/S0006-3495(61)86902-6.
13. Terman D. An Introduction to Dynamical Dystems and Neuronal Dynamics. In: Tutorials in Mathematical Biosciences I. Lecture Notes in Mathematics. Berlin, Heidelberg: Springer-Verlag; 2005. P. 21-68.DOI: 10.1007/978-3-540-31544-5_2.
14. Glyzin SD, Kolesov AYu. On a Method of Mathematical Modeling of Electrical Synapses. Diff Equat. 2022;58:853-868. DOI: 10.1134/S0012266122070011.
15. Kolesov AYu, Mishchenko EF, Rozov NKh. Relay with delay and its С1-approximation. Proc. Steklov Inst. Math. 1997;216:119-146.
16. Glyzin SD, Kolesov AYu. Traveling waves in fully coupled networks of linear oscillators. Comput. Math. Math. Phys. 2022;62(1):66-83. DOI: 10.1134/S0965542522010079.
17. Glyzin SD, Kolesov AYu, Rozov NKh. Buffering in cyclic gene networks. Theoret. and Math. Phys. 2016;187(3)935-951. DOI: 10.1134/S0040577916060106.
18. Brown PN, Byrne GD, Hindmarsh AC. VODE: A Variable Coefficient ODE Solver // SIAM J. Sci. Stat. Comput. 1989;10(5):1038-1051.DOI: 10.1137/0910062.
19. Glyzin SD, Kolesov AYu. Periodic two-cluster synchronization modes in fully coupled networks of nonlinear oscillators. Theor. Math. Phys. 2022;212:1073-1091. DOI: 10.1134/S0040577922080049.
Глызин Дмитрий Сергеевич — родился в Ярославле (1983), окончил Ярославский государственный университет (2004). Защитил диссертацию на соискание ученой степени кандидата физико-математических наук (2006). Область научных интересов: нелинейная динамика, математическое моделирование.
Россия, 150003 Ярославль, ул. Советская, 14
Ярославский государственный университет имени П. Г. Демидова E-mail: [email protected] ORCID: 0000-0002-0701-622X AuthorlD (eLibrary.Ru): 749418
Глызин Сергей Дмитриевич — родился в Новгородской области (1957), окончил Ярославский государственный университет (1979). Защитил диссертации на соискание ученой степени кандидата физико-математических наук в ВГУ (1991) и доктора физико-математических наук в НИЯУ МИФИ (2009) в области теории нелинейных колебаний. Профессор, заведующий кафедрой компьютерных сетей. Опубликовал более 100 научных работ. Приказом Министерства образования и науки награжден нагрудным знаком «Почетный работник высшего профессионального образования». Заместитель главного редактора научного журнала «Моделирование и анализ информационных систем».
Россия, 150003 Ярославль, ул. Советская, 14
Ярославский государственный университет имени П. Г. Демидова E-mail: [email protected], [email protected] ORCID: 0000-0002-6403-4061 AuthorID (eLibrary.Ru): 35140
Колесов Андрей Юрьевич — родился в городе Старый Оскол Белгородской области (1961). Окончил Ярославский государственный университет (1984). Защитил диссертации на соискание ученой степени кандидата физико-математических наук (1988) и доктора физико-математических наук (1991) в Математическом институте имени В. А. Стеклова РАН в области дифференциальных уравнений и динамических систем. Профессор кафедры дифференциальных уравнений. Опубликовал более 200 научных работ, среди которых 6 монографий.
Россия, 150003 Ярославль, ул. Советская, 14
Ярославский государственный университет имени П. Г. Демидова E-mail: [email protected] ORCID: 0000-0001-5066-0881 AuthorlD (eLibrary.Ru): 5544
