CC BY
5
УДК 530.182
ДИНАМИКА ЦЕПОЧКИ ВЗАИМОДЕЙСТВУЮЩИХ
НЕЙРОНОВ С НЕЛОКАЛЬНОЙ СВЯЗЬЮ, ЗАДАННОЙ
ОПЕРАТОРОМ ЛАПЛАСА ДРОБНОГО И ПЕРЕМЕННОГО
ПОРЯДКОВ С НЕЛИНЕЙНЫМИ ФУНКЦИЯМИ МОДЕЛИ
ХИНДМАРШ-РОУЗ
И. С. Фатеев, А. А. Полежаев
Предложена модель, описывающая динамику одномерной цепочки взаимодействующих нейронов, объединенных нелокальной связью, основанная на свойствах дробного лапласиана, как с постоянным, так и с переменным показателем. В результате численного моделирования распространения потенциала действия для различных типов нелокальных взаимодействий (с постоянным набором параметров, соответствующих нелинейным функциям модели Хиндмарш-Роуз, а также коэффициентов диффузии) установлено возникновение различных пространственно-временных режимов. При заданных параметрах, в случае классической диффузии наблюдается общая тенденция к синхронному поведению. В супердиффузионном случае формируются режимы кластерного возбуждения потенциала действия. Для систем с переменным показателем дробного лапласиана характерна пространственная анизотропия возникающих структур. Продемонстрировано, что внедрение дально-действующих связей, реализуемых в системе введением дробного оператора Лапласа, способно предоставить дополнительные возможности описания динамических свойств взаимодействующих нейронов.
Ключевые слова: коллективная динамика нейронов, модель Хиндмарш-Роуз, дробный оператор Лапласа, дифференциальные операторы с переменным порядком. ФИАН, 119991 Россия, Москва, Ленинский пр-т, 53; e-mail: i.fateev@lebedev.ru, apol@lpi.ru.
Введение. Совокупность взаимодействующих нейронов представляет собой сложную динамическую систему с богатыми внутренними свойствами. Разнообразные характеристики таких систем, начиная от частотных параметров активации конкретных нервных клеток, заканчивая пространственно-временной динамикой их совокупностей, описываются системами нелинейных дифференциальных уравнений. В частности, для описания распределенных систем используют реакционно-диффузионные уравнения и соответствующее им шаблоны разностной аппроксимации (см., напр., [1-4]). Было показано, что в рамках систем с данной кинетикой могут возникать спирали и гексагональные структуры, которые объясняют некоторые особенности коллективной динамики локально (элемент взаимодействует лишь с соседними элементами) взаимодействующих нервных клеток [5, 6]. Обоснование диффузионного подхода, в такой постановке, отождествляют с переносом ионов К + [7]. Не менее важным инструментом описания систем нелокально (элемент взаимодействует с несколькими другими элементами) или глобально (элемент взаимодействует со всеми элементами) взаимодействующих нейронов являются интегро-дифференциальные уравнения [8-10]. Конкретные варианты их применения можно найти в работах [11, 12]. Альтернативой представленным континуальным подходам может служить множество дискретных систем, в которых роль диффузии и интегральных операторов отдается операторам суммирования в сочетании с определенным типом взаимосвязи между элементами [13-16].
Данные методы уже принесли многочисленные интересные результаты и для одномерных систем. Так, в цепочках взаимодействующих нейронов были обнаружены состояния, называемые химерами, сочетающие в себе как области с высокой степенью когерентности, так и инкогерентные участки [17-20]. Важной особенностью таких состояний является пространственная неоднородность возникающих структур. Эмпирические данные, которые соответствуют химерным состояниям, представлены в [21]. Аналогично многомерным моделям, одномерные системы представляют особый интерес и дают качественное понимание процессов коллективной нейронной активности [22, 23].
Последние исследования демонстрируют наличие связи между групповой активностью нейронов и супердиффузионным кинетическим механизмом [24, 25], описание которого производится в терминах полетов Леви. В таком ключе, особую роль приобретают степенные законы, которые также связывают с универсальным механизмом образования нейронных лавин [26], вероятностно топологическими характеристиками развитых крупномасштабных сетей [27], а также с некоторыми закономерностями в социально-поведенческих приложениях [28]. Эти данные могут служить предпосылка-
ми для применения иного подхода к изучению кинетики распространения потенциала действия в нервной ткани, основанного на дробно-дифференциальном подходе.
В рамках данного исследования мы рассматриваем пространственно-временную динамику в одномерном случае для совокупности взаимодействующих нервных клеток. Для описания нелокальных связей в качестве основы модели мы выбираем аппарат супердиффузии с дробным оператором Лапласа, постоянного и переменного порядков, а также сравниваем данную систему с классической, кинетика которой задается стандартным лапласианом и его трехточечной схемой аппроксимации, формирующей собой локальный тип соединения. Ввиду общей математической лаконичности и богатого внутреннего функционала, способного описать практически все известные примеры динамического поведения точечного нейрона [29], мы используем нелинейные функции, соответствующие модели Хиндмарш-Роуз. Для сравнения различных типов соединений (локального и нелокального) был использован единый набор параметров, соответствующий однородному режиму пульсаций классической системы.
Цель текущего исследования в демонстрации возможностей дробно-дифференциального подхода с постоянным и переменным порядком пространственной производной для описания коллективной динамики нелокально взаимодействующих нервных клеток в одномерном случае. Мы показываем, что показатель дробного лапласиана может являться дополнительным параметром, влияющим на коллективную динамику нелокально связанных нейронов. Более того, его переменная форма определяет топологию соединений и может использоваться в качестве универсального механизма организации синаптических связей.
Одномерная модель супердиффузионной связи нейронов.
А. Одномерная модель с постоянным порядком. Совокупность взаимодействующих нервных клеток может быть описана системой нелинейных дифференциальных уравнений, модифицированных путем добавления некоторого оператора, ответственного за тип связи:
Конкретные реализации данного оператора, помимо уже упомянутых работ [8-10], можно обнаружить также в [30]. В свою очередь, во многих задачах современной физики и биологии для описания нелокального поведения исследуемых величин часто используется дробно-дифференциальный подход [31]. В данном случае, нелокальность может быть обусловлена степенной асимптотикой в поведении исследуемых решений. Для опи-
(1)
сания подобных систем пользуются обобщением классического оператора Лапласа на его дробный аналог [32, 33].
/ ' V . (0,1), (2)
К,
а
где ' - размерность пространства, V - показатель дробного лапласиана, Г(х) - гамма-функция. Данный оператор является частным случаем (1) с заданным степенным ядром. Иное представление оператора (2) может быть записано в терминах преобразования Фурье следующим образом [31-33]:
г|(-А)а/2/(г")|(к) = \к\а/(к). (3)
Прямое указание на супердиффузионный кинетический механизм [24, 25], связь степенных функций с нейронными лавинами [26] и вероятностно-топологическими особенностями крупных сетей [27], а также следствие дробно-дифференциального подхода из общей концепции нелокального взаимодействия, свидетельствуют о потенциале моделей с дробным оператором Лапласа для задач нейронаук. Более того, данный метод способен предоставить дополнительную связь между нелокальным и глобальным подходами, основанными на применении интегро-дифференциальных уравнений и реакционно-диффузионными системами, ответственными, как правило, за локальный тип взаимодействия.
На основе данных определений, для описания нелокальных свойств динамических сетей взаимодействующих нейронов, мы вводим следующую одномерную трехкомпо-нентную систему:
'п(х,1)/(И = —Би(—Ах)аи/2п(х,1) + у(х,1) — ап3(х,Ь) + Ьп2(х,Ь) — т(х,Ь) + 1еХ ¿у(х,г)/'г = —Б1] (—Ах)а"/2у(х,г) + с — (п2 (х,г) — у(х,г) (т(х, Ь)/(И = —Бт(—АхХ)ат/2т(х, Ь) + т^в(п(х, Ь) — п0) — т(х, ,
(4)
в которой п(х,Ь) - потенциал действия, у(х,Ь) и т(х,Ь) описывают динамику переноса ионов К + и Са+ через внешнюю мембрану нервной клетки в точках х в момент времени Ь, 1ехЬ - внешний стимул, Пи,Пь,От - коэффициенты диффузии соответствующих компонент (Ои ^ ^; Ди ^ От). Наличие дробной диффузии во всех уравнениях системы (4) является обобщением идеи обыкновенной диффузии компонент V и т при
наличии супердиффузионного режима для компоненты и. Аналогичный подход развивается, в частности, в работах [34-38]. Нелинейные функции, стоящие в правых частях, соответствуют модели Хиндмарш-Роуз и качественно описывают локальную динамику отдельного нейрона. Одномерный оператор Лапласа, в предположении об отсутствии влияния областей, находящихся за пределами рассматриваемого диапазона [а, Ь], может быть представлен в следующей форме:
1 " "
(-Ах)а/21 (х)
2соъ(па/2)
(х) +х (х)
(5)
где правосторонняя „В^ и левосторонняя хВа дробная производная задаются соотношениями:
х
в» !(х) —_1_дИ+1 /" !(х^(х' и> 0 (6)
аДх/ (х) = Г(1 -{и})Вх у (х - х/)М , и> 0 (6)
а
Ь
и» г(х) — (-1)м+1 вм+1 Г !(х')Лх1 и> 0 (7)
хВь /(х)-г(1 -{V})Вх ] (х/ - х)М , и> 0- (7)
х
Тут {V} - дробная часть показателя, а [и] - целая [32, 33]. При показателе а — 2 мы переходим к производной второго порядка и, как следствие, к системе классических нелинейных реакционно-диффузионных уравнений.
В. Одномерная модель с переменным порядком. Дробно-дифференциальные уравнения с переменным порядком производных являются актуальным инструментом описания многочисленных явлений в задачах физики и биологии [39]. На данный момент уже существуют исследования, которые выделяют конкретные преимущества данного подхода [40]. Мы предлагаем использовать аппарат дробного дифференцирования переменного порядка для выделения пространственных областей с различными типами межклеточного взаимодействия. Система дифференциальных уравнений в таком случае может быть представлена:
(и(х,г)/(г — -в„(-Дх)аи(х,4)/2и(х,г) + у(х,г) - аи3(х,г) + Ьи2(х,г) - ш(х,г) + 1ехь (у(х,г)/М — -В,и(-Дх)^(х'*)/2у(х,г) + с - (!и2(х,г) - у(х,г) (ш(х, Ь)/(Ь — -Пт(-Дх)ат(х'*)/2ш(х, Ь) + г^в(и(х, Ь) - и0) - ш(х, Ь)^ ,
(8)
где лапласиан переменного порядка задается аналогично (3) в терминах преобразования Фурье [41]. Предполагается, что данный подход способен предоставить дополнительные возможности для описания сложно организованных сетей взаимосвязанных нейронов.
Модель (8), в случае постоянных коэффициентов диффузии, сводится к явной итерационной схеме, представленной в [42]:
, ,1Ч г+1 ,Л /гл Хтах
»+1 , Л
п^1 = п* + БиТ^ Е д^,*<-+1 + ПиТ% Е п3 + '1/1 (п*, т*)
3=0 3=0
»+1 ,Л Хтах
»+1 /л
vЧ+1 = V* + Б,г(1 Е д^гч3 + БгЩ Е 9(-)1,чVI-, + /(п*, V*, т*)
3=0 3=0
»+1 ,л /^.ч Хтах
»+1 , л
т*+1 = т* + Д„г$ Е 9(+)1, * ^-,+1 + А^З Е 9 й,* т3 + Л/э(п* , т*).
3=0 3=0
(9)
1.0 0.5 0.0 -0.5
-1.0
-1.5
-2.0
_0 2 4 6 8 И)
Рис. 1: Основная часть: весовые функции д3 = д3 (для различных постоянных показателей а дробного лапласиана (—А)а/2), фигурирующие в (9). Вставка: схематичное представление воздействия одномерного дробного оператора Лапласа и его связь с особенностями соединения для 1-го элемента динамической системы. Линии, отображенные пунктиром, соответствуют локальному соединению (аЧ-1 = аЧ+1 = 2). Сплошные линии характеризуют супердиффузионный режим с аЧ-1 = аЧ+1 < 2.
Функции /1 (и, V, т), /2(и, V, т), /3(и, V, т) соответствуют нелинейным функциям модели Хиндмарш-Роуз. Множители перед оператором суммы определены как г(1к =
sec(nak/2)
dt(dx)
-ak+i и r(2) i,k
sec(nak/2) ,_nk
—dt(dx) "i-1. Второе и третье слагаемые в
2 4 ' 2 правой части есть конечные разности для левосторонней и правосторонней производной, соответственно из (5). Весовые функции определяются следующим соотношением:
g(j) gi,k
(-1)J
aM = (-I)j r(ak + 1) j r(ak - j + 1)r(j + 1)'
:iq)
Рассчитанные д^ для постоянных значений а, а также схематичное представление воздействия одномерного дробного лапласиана и его влияние на взаимосвязь элементов системы представлены на рис. 1. В случае, когда ак = 2 переходим к классической явной схеме решения диффузионного уравнения.
Рис. 2: (а) Поведение точечной системы при заданных параметрах с начальными условиями Мшй = и* — 1, Vinit = V* + 1, ттй = т* — 10. Ось абсцисс соответствует времени Ь, а ось ординат - компонентам и^,т. (Ь) Динамика точечной системы в пространстве компонент и^,т со сценарием перехода к стационарному режиму.
Результаты. Для проведения численных расчетов были использованы следующие значения параметров системы (4): а =1, b = 3, c =1, d = 5, r = 0.008, s = 4, u0 = -1.6, lext = 1.315, при которых точечная система имеет устойчивую стационарную точку. На рис. 2 показана динамика возврата системы к стационару после начального отклонения. Фрагмент (а) демонстрирует поведение компонент u,v,m точечного нейрона при заданных параметрах нелинейной модели. На фрагменте (b) представлена динамика то-
Рис. 3: 1) расчетные данные для классической модели (а» = 2). 2), 3) расчетные данные для модели (4) с показателями аи = 1.8 и аи = 1.2, соответственно, при классической схеме распространения двух других компонент (а, = 2,ат = 2). 4) описывает, пространственно-временную динамику в случае аи = 1.5 и а, = 1.8, ат = 1.8. Область
(a)зона, в которой имеется синхронизация, при внедрении нелокальных соединений;
(b) зона с нарушением синхронизации; (с) зона перехода системы к стационару. Горизонтальные пунктирные линии внедрены для удобства сопоставления активационного процесса.
чечной системы в фазовом пространстве и ее сходимость к устойчивому стационарному состоянию {п*^*,т*}.
Моделирование распределенной системы было организовано на последовательности из 100 точек с пространственно-временной дискретизацией: 'х = 0.005, (Ь = 0.05, а также разделением по физическим процессам, где кинетическая часть рассчитывалась
по явной итерационной схеме (9), а нелинейная обратная связь методом Рунге-Кутты 4-го порядка. Также, были выбраны следующие значения коэффициентов диффузии: Ои = 10-4, = 10-5, От = 10-5. В качестве начальных условий во всей пространственной области были заданы значения, соответствующие стационарному состоянию точечной системы: и* = —1.317, V* = —7.678, т* = 1.130. Несколько центральных элементов были выведены из стационарного состояния: и^ = и* — 1, Vinit = V* + 1, т[пц = т* — 10. Для исключения возникновения вторичного динамического процесса с образованием дополнительных волн, а также для демонстрации развития динамики исключительно центральной части, в качестве краевых условий были использованы стационарные значения.
Расчетные данные для различных комбинаций аи,а^,ат представлены на рис. 3. Фрагмент 1): При заданных параметрах, в рамках классической локальной схемы диффузионной связи (аи = а.и = ат = 2) наблюдается режим общей синхронизации процесса активации нейронов ввиду повышенной связности локально взаимодействующих элементов. Фрагмент 2): Внедрение нелокальных связей (аи = 1.8, а^ = ат = 2) приводит к нарушению синхронизации и, как следствие, к иным пространственно-временным структурам, по сравнению с классической схемой взаимодействия (см. области 1, Ь) и 2, Ь)). Фрагмент 3): Уменьшение показателя аи (аи = 1.2, а^ = ат = 2) приводит к более ранней десинхронизации системы (см. область 3, а)). Фрагмент 4): Супердиффузионный режим для всех трех компонент (аи = 1.5, а^ = ат = 1.8) приводит к ещё большей рассинхронизации и формированию кластерной структуры. Рис. 3 демонстрирует общую тенденцию запаздывания перехода системы к стационарному состоянию (см. область с) для всех изображений.
Вычисления, при заданных параметрах нелинейной модели, а также коэффициентах диффузии, проделаны для конкретных конфигураций связей взаимодействующих нейронов и значений показателей аДж), соответственно. Тем самым, система была разделена на области с различным типом соединений между элементами. Данные представлены на рис. 4. Полученные результаты демонстрируют общие черты анизотропного динамического поведения системы при ее переходе к стационарному состоянию (фрагмент 1), область а)). При введении переменной связи по всем компонентам повышенная кластеризация проявляется анизотропно (фрагмент 2), область Ь)). Наконец, данные, представленные на рис. 5, демонстрируют пространственно-временную динамику при изменении типа соединений между элементами, в случае, когда дальнодействующие связи постепенно заменяются на локальные. Данное изменение приводит к увеличению
Рис. 4: Расчетные данные для модели (9) с конкретными формами зависимостей а»(х), представленными справа от диаграммы. Параметры 1) рисунка: аи = аи(х), а, = ат = 2. Для 2) характерны зависимости а»(х) для всех трех компонент.
масштаба структур (см. область внутри пунктирной линии). Сам процесс может быть отождествлен с изменением синаптических связей между нейронами (подробнее в [43]).
Характерным общим отличием супер диффузионной модели, как и модели с переменным дробным показателем, является степенной спад компонент п, V, т относительно возмущенного состояния. В дополнении отметим, что повышенная кластеризация в супердиффузионном случае и в случае систем с переменным порядком, может быть связана с масштабными преобразованиями структур, которые наблюдаются в системах с дробными производными (см., напр., [44-46]). Мы полагаем, что совокупность параметров, ответственных за нелокальность (а») и силу взаимной связи между элементами (Б»), а также их взаимная организация, могут служить адекватными инструментами
Рис. 5: Расчетные данные для модели (9) с динамической зависимостью аи(х,^) при заданных а,0 = ат = 2. Область внутри пунктирной .линии демонстрирует изменение масштаба структур, вызванное изменением типа соединений.
для формирования сложных топологий сетей, основанных на универсальном законе взаимодействия (2).
зЗаключение. В данной работе была рассмотрена пространственно-временная динамика связанных нейронов с учетом различного характера промежуточных нелокальных взаимодействий между ними, организованная с использованием дробного оператора Лапласа постоянного и переменного порядков. Произведено сравнение классического диффузионного переноса (а = 2) с супердиффузионным (а < 2), а также разобран случай с функциональной зависимостью аДх) и аи(х,Ь) при равных общих параметрах, соответствующих модели Хиндмарш-Роуз и коэффициентах диффузии Показано, что при единых начальных и граничных условиях система, с внедренным су-
пердиффузионным взаимодействием, демонстрирует иное пространственно-временное поведение на больших временах. Для совокупности параметров, при которых классическая система (а» = 2) демонстрирует режим общей синхронизации, системы с постоянными показателями а» < 2 воспроизводят переход к кластерному режиму возбуждения, а пульсации поддерживают больший промежуток времени. Если кинетика задается несколькими супердиффузионными уравнениями, темпы развития кластерных структур выше. При конкретных формах показателей а(х) проявляется анизотропия пространственно-временных структур в рассматриваемом масштабе. Системы с пространственно-временной зависимостью показателя аи(х,Ь) способны описать поведение изменчивой сети взаимодействующих элементов. В тех областях пространства, в которых задано нелокальное взаимодействие, происходит развитие структур на больших временах.
На наш взгляд, параметры, ответственные за нелокальность взаимодействия между элементами системы в общем и супердиффузионный подход в частности, способны предоставить дополнительные возможности для описания режимов возбуждения цепочки связанных нейронов. В работе показано, что вышеизложенное справедливо даже для тех совокупностей параметров нелинейных функций, для которых в случае обычной диффузионной связи имеет место синхронное поведение.
Работа выполнена при финансовой поддержке РНФ (грант № 22-21-00546).
ЛИТЕРАТУРА
[1] S. K. Sharma, A. Mondai, A. Mondai, et al., International Journal of Bifurcation and Chaos 31, 1 (2021). DOI: 10.1142/S0218127421300305.
[2] A. Mondai, R. K. Upadhyay, A. Mondai, S. K. Sharma, Applied Mathematics and Computation 423, 1 (2022). DOI: 10.1016/j.amc.2022.127010.
[3] K. Wang, Z. Teng, H. Jiang, Communications in Nonlinear Science and Numerical Simulation 17, 3866 (2012). DOI: 10.1016/j.cnsns.2012.02.020.
[4] M. Jun, Y. He-Ping, L. Yong, L. Shi-Rong, Chinese Physics B 18, 98 (2009). DOI: 10.1088/1674-1056/18/1/017.
[5] A. Keane, P. Gong, The Journal of Neuroscience 35, 1591 (2015). DOI: 10.1523/JNEUROSCI.1669-14.2015.
[6] R. G. Townsend, P. Gong, PLOS Computational Biology 14, 1 (2015). DOI: 10.1371/journal.pcbi.1006643.
[7] D. M. Durand, E. Park, A. L. Jensen, Phil. Trans. R. Soc. B 365, 2347 (2010). DOI: 10.1098/rstb.2010.0050.
[8] Si. Amari, Biological cybernetics 27, 77 (1977). DOI: 10.1007/BF00337259.
[9] S. I. Amari, Bulletin of mathematical biology 42, 339 (1980). DOI: 10.1016/S0092-8240(80)80055-3.
[10] B. Ermentrout, Reports on progress in physics 61, 353 (1998). DOI: 10.1088/00344885/61/4/002.
[11] K. A. Richardson, S. J. Schiff, B. J. Gluckman, Physical review letters 94, 028103 (2005). DOI: 10.1103/PhysRevLett.94.028103.
[12] X. Huang, W. C. Troy, Q. Yang, et al., Journal of Neuroscience 24, 9897 (2004). DOI: 10.1523/JNEUROSCI.2705-04.2004.
[13] Z. T. Njitacke, S. S. Muni, S. Seth, et al., Physica Scripta 98, 045210 (2023). DOI: 10.1088/1402-4896/acbdd1.
[14] S. S. Muni, Z. Njitacke, C. Feudjio, et al., Chaos Theory and Applications 4, 119 (2022). DOI: 10.51537/chaos.1144123.
[15] D. Linaro, M. Righero, M. Biey, M. Storace, Synchronization properties in networks of Hindmarsh-Rose neurons and their PWL approximations with linear symmetric coupling. 2009 IEEE International Symposium on Circuits and Systems, pp. 1685-1688. DOI: 10.1109/ISCAS.2009.5118098.
[16] Q. Zheng, J. Shen, Applied Mathematics and Computation 381, 125304 (2020). DOI: 10.1016/j.amc.2020.125304.
[17] J. Hizanidis, V. G. Kanas, A. Bezerianos, T. Bountis, International Journal of Bifurcation and Chaos 24, 1 (2014). DOI: 10.1142/S0218127414500308.
[18] B. K. Bera, D. Ghosh, M. Lakshmanan, Physical Review E 93, 012205 (2016). DOI: 10.1103/PhysRevE.93.012205.
[19] Z. Wei, F. Parastesh, H. Azarnoush, et al., Europhysics Letters 123, 48003 (2018). DOI: 10.1209/0295-5075/123/48003.
[20] Z. Wang, Y. Xu, Y. Li, et al., Chaos, Solitons and Fractals 148, 110976 (2021). DOI: 10.1016/j.chaos.2021.110976.
[21] S. Majhi, B. K. Bera, D. Ghosh, Physics of life reviews 28, 100 (2019). DOI: 10.1016/j.plrev.2018.09.003.
[22] A. Mondal, S. K. Sharma, R. K. Upadhyay, et al., Phys. Rev. E 99, 1 (2019). DOI: 10.1103/PhysRevE.99.042307.
[23] A. Mondal, A. Mondal, S. K. Sharma, R. K. Upadhyay, Nonlinear Dynamics 105, 3515 (2021). DOI: 10.1007/s11071-021-06787-4.
[24] G. Chen, P. Gong, Science Advances 8, 1 (2022). DOI: 10.1126/sciadv.abl4995.
[25] Y. Qi, P. Gong, Nature communications 13, 1 (2022). DOI: 10.1038/s41467-022-32279-z.
[26] J. M. Beggs, D. Plenz, The Journal of Neuroscience 23, 11167 (2003). DOI: 10.1523/JNEUR0SCI.23-35-11167.2003.
[27] A. L. Barabasi, R. Albert, Science 286, 509 (1999). DOI: 10.1126/science.286.5439.509.
[28] A. Baronchelli, F. Radicchi, Chaos, Solitons and Fractals 56, 101 (2013). DOI: 10.1016/j.chaos.2013.07.013.
[29] M. Storace, D. Linaro, E. de Lange, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 033128 (2008). DOI: 10.1063/1.2975967.
[30] F. Parastesh, S. Jafari, H. Azarnoush, et al., Physics Reports 898, 1 (2021). DOI: 10.1016/j.physrep.2020.10.003.
[31] R. Klages, G. Radons, I. M. Sokolov, Anomalous Transport. Foundations and Applications (WILEY-VCH Verlag GmbH and Co. KGaA, 2008).
[32] В. В. Учайкин, Метод дробных производных (Ульяновск, Артишок, 2008).
[33] S. G. Samko, A. A. Kilbas, O. I. Marichev, Fractional Integrals and Derivatives. Theory and Applications (Gordon and Breach, Yverdon, 1993).
[34] Y. Wu, J. Xu, D. He, D. J. Earn, Chaos, Solitons and Fractals 23, 1605 (2005). DOI: 10.1016/j.chaos.2004.06.077.
[35] X. Li, C. Han, Y. Wang, Fractal and Fractional 6, 136 (2022). DOI: 10.3390/fractalfract6030136.
[36] F. Liu, I. Turner, V. Anh, et al., Anziam Journal 54, 608 (2012). DOI: 10.21914/anziamj.v54i0.6372.
[37] F. Liu, P. Zhuang, I. Turner, et al., Journal of Computational Physics 293, 252 (2015). DOI: 10.1016/j.jcp.2014.06.001.
[38] Q. Zheng, J. Shen, Applied Mathematics and Computation 381, 1 (2020). DOI: 10.1016/j.amc.2020.125304.
[39] S. Patnaik, J. P. Hollkamp, F. Semperlotti, Proc. R. Soc. A 476, 1 (2020). DOI: 10.1098/rspa.2019.0498.
[40] H. G. Sun, W. Chen, H. Wei, Y. Q. Chen, The european physical journal special topics 193, 185 (2011). DOI: 10.1140/epjst/e2011-01390-6.
[41] M. D. Ruiz-Medina, V. V. Anh, J. M. Angulo, Stochastic Analysis and Applications 22, 775 (2004). DOI: 10.1081/SAP-120030456.
[42] P. Zhuang, F. Liu, V. Anh, I. Turner, Siam J. Numer. Anal. 3, 1760 (2009). DOI: 10.1137/080730597.
[43] G. Mongillo, S. Rumpel, Y. Loewenstein, Current opinion in neurobiology 46, 7 (2017). DOI: 10.1016/j.conb.2017.06.006.
[44] H. G. Lee, Journal of Computational and Applied Mathematics 333, 395 (2018). DOI: 10.1016/j.cam.2017.09.007.
[45] A. Bueno-Orovio, K. Burrage, Communications in Nonlinear Science and Numerical Simulation, 107120 (2023). DOI: 10.1016/j.cnsns.2023.107120.
[46] K. M. Owolabi, B. Karaagac, D. Baleanu, Mathematical Methods in the Applied Sciences 44(5), 4018 (2021). DOI: 10.1002/mma.7007.
Поступила в редакцию 13 февраля 2023 г.
После доработки 4 мая 2023 г. Принята к публикации 5 мая 2023 г.
