Численный анализ стационарных решений систем с запаздывающим аргументом в математической иммунологии

Полный текст

1. Введение Математическое моделирование является важнейшим аналитическим инструментом современной математической иммунологии [7, 13]. Интенсивные исследования в области инфекции вирусами SARS-CoV-2 последних лет демонстрируют возможности математических моделей и возникающие задачи по развитию вычислительных методов и технологий их исследования [10]. Фундаментальной по значимости проблемой современной иммунологии является изучение механизмов развития неблагоприятных форм течения вирусных инфекций, в частности, хронического вирусного гепатита В. Известно, что хронический вирусный гепатит является многофакторным заболеванием [8]. Ранее нами путем анализа чувствительности решений математической модели противовирусного иммунного ответа было установлено, что изменение в параметрах активации антигенпрезентирующей функции макрофагов может служить одним из факторов перехода от острой формы инфекции с выздоровлением, к хронической персистенции вирусов [4]. Вместе с тем, систематическое исследование характеристик хронического вирусного гепатита В остается недостаточно изученной проблемой. В частности, необходимо находить области в пространстве параметров математических моделей, отвечающие бистабильным или мультистабильным режимам взаимодействия инфекции и иммунной системы организма, а также изучать возможности переходов между ними. Ключевым для анализа данных свойств является наличие эффективных методов вычисления и численного анализа устойчивости стационарных решений нелинейных систем дифференциальных уравнений с запаздывающим аргументом, используемых для моделирования вирусных инфекций. Бистабильность, как возможность системы «вирус - организм человека» сосуществовать в двух устойчивых равновесных состояниях, является чрезвычайно важным свойством, поскольку позволяет решать задачу функционального лечения вирусной инфекции путем перехода из хронического состояния с более высокой вирусной нагрузкой в более благоприятное состояние с низкой вирусной нагрузкой за счет активации компонент иммунной системы. Биологические аргументы в подтверждение реалистичности данной концепции представлены в работе [3]. Следует отметить, что перевод инфекционного процесса в состояние с более низкой вирусной нагрузкой является нетривиальной практической задачей, для решения которой необходима активация двух звеньев иммунитета - гуморального и Т-клеточного. При этом мультипликативный эффект вклада антител и цитотоксических Т лимфоцитов в снижение уровня инфекционного процесса позволяет снизить пороги численности данных популяций, требуемые для контроля инфекции на более низком уровне. Изучение принципиальных условий бии (или) мультистабильности инфекционного процесса и кооперативности звеньев иммунитета в обеспечении контроля уровня вирусной нагрузки возможно провести на основе математической модели вирусного заболевания с применением методов бифуркационного анализа. Математическая модель противовирусного иммунного ответа, разработанная Г.И. Марчуком и Р.В. Петровым [14], представляет собой эффективный инструмент для исследования условий сосуществования стационарных состояний хронической вирусной инфекции, поскольку обладает необходимой для этого полнотой описания реакций адаптивного иммунитета и калибрована по данным острой инфекции вирусами гепатита В, см. [15]. Ранее с помощью данной модели были изучены отдельные кинетические механизмы хронизации вирусного гепатита В по данным двух вариантов динамики инфекции (острой и хронической) у добровольцев, см. [4]. На основе модели был предсказан мультипликативный эффект антител и цитотоксических Т-клеток в обеспечении защиты организма от вирусной инфекции для гриппа А, вирусного гепатита В и валидирован по данным инфекции ВИЧ, см. [2]. Сформулированная в виде среднеразмерной системы дифференциальных уравнений с запаздывающим аргументом, данная модель является полилинейной функцией от n аргументов (n :( 4). Это предполагает возможность существования до четырех нетривиальных положений равновесия и, следовательно, наличие мультистабильности в пространстве состояний модели. Предварительное исследование существования стационарных решений в модели Марчука-Петрова было выполнено нами ранее в работах [12, 18]. Целью настоящей работы является развитие методов численного поиска корней систем нелинейных алгебраических уравнений, соответствующих положениям равновесия модели, и анализа их локальной устойчивости. С помощью разработанной технологии будет проведено полноценное исчерпывающее исследование условий сосуществования нескольких стационарных решений и переходов между ними, что необходимо для формирования механизменного понимания различных режимов стационарной динамики хронического вирусного гепатита В. Работа состоит из шести разделов. В разделе 2 кратко описана рассматриваемая математическая модель противовирусного иммунного ответа Марчука-Петрова. В разделе 3 описан и обоснован предлагаемый метод вычисления и трассирования стационарных решений по параметрам модели. Результаты вычисления стационарных решений модели Марчука-Петрова и анализа зависимости вычисленных решений от параметров модели приводятся в разделе 4, где также дается их содержательная интерпретация. Итог работы подводится в заключительном разделе 5. 2. Математическая модель противовирусного иммунного ответа Марчука-Петрова Математическая модель противовирусного иммунного ответа Марчука-Петрова сформулирована в [4, 14] в виде системы нелинейных дифференциальных уравнений с запаздывающим аргументом. Система описывает скорость изменения во времени концентрации следующих популяций: вирусных частиц Vf ; зараженных вирусами клеток органа-мишени CV ; разрушенных клеток органа-мишени m; антигенпрезентирующих клеток (макрофагов) MV ; CD4+ T-лимфоцитов - помощников клеточного иммунитета (Th1) HE ; CD4+ Т-лимфоцитов - помощников гуморального 688 М. Ю. ХРИСТИЧЕНКО и др. иммунитета (Th2) HB ; CD8+ Т-лимфоцитов - киллеров E, уничтожающих зараженные вирусами клетки; В-лимфоцитов B; плазматических клеток P, производящих антитела; антител F, нейтрализующих вирусы. Будем предполагать далее, что переменные Vf и F имеют размерность частиц/мл, а остальные переменные - клеток/мл. Уравнения модели построены на основе балансных соотношений для скоростей процессов рождения, дифференцировки и гибели компонентов системы «вирус - организм хозяина». Первые три уравнения модели описывают кинетику распространения инфекции в клетках органа-мишени с учетом действия факторов иммунной защиты. В уравнении для вирусной популяции, наряду с секрецией инфекционных вирусных частиц зараженными клетками в течение их функционирования (первое слагаемое в первом уравнении), описывается процесс высвобождения вирусных частиц (например, в случае вирусного гепатита В речь идет о 22-нм частицах HBsAg) в результате разрушения инфицированных клеток CD8+ T-лимфоцитами-киллерами. С помощью функции ξ(m) учитывается ослабление иммунного ответа вследствие разрушения зараженных клеток органа-мишени. Активация CD4+ T-лимфоцитов-помощников происходит в результате взаимодействия с антигенпрезентирующими клетками MV , а активация CD8+ T-лимфоцитов и В-лимфоцитов предполагает наличие сигналов от антигенпрезентирующих клеток и CD4+ Tлимфоцитов-помощников. Соответственно, скорости этих процессов описывается билинейными и трилинейными функциями с запаздывающим аргументом, учитывающим длительность прохождения делящимися клетками всех фаз клеточного цикла. Остальные взаимодействия параметризованы с применением закона действующих масс, схемы «хищник - жертва» и экспоненциальной кинетики естественного убывания численности антител, вирусов и клеток. Система уравнений модели имеет следующий вид: d dt Vf (t) = νCV (t)+ nbCECV (t)E(t) - γVF Vf (t)F (t) - γVM Vf (t) - - γVC Vf (t) fC0 - CV (t) - m(t)l , d dt CV (t) = σVf (t) fC0 - CV (t) - m(t)l - bCECV (t)E(t) - bmCV (t), d d dt m(t) = bCECV (t)E(t)+ bmCV (t) - αmm(t), dt MV (t)= γMV M 0Vf (t) - αM MV (t), d HE (t) = bE [ξ(m)ρE MV (t - τ E )HE (t - τ E ) - MV (t)HE (t)] - dt H H H H - bHE MV (t)HE (t)E(t)+ αE (H0 - HE (t)), p H E d (2.1) p dt E(t) = bE [ξ(m)ρEMV (t - τE )HE (t - τE )E(t - τE ) - MV (t)HE (t)E(t)] - § bECCV (t)E(t)+ αE (E0 - E(t)), bB [ξ(m)ρB MV ( H H t - τ B )HB (t - τ H B ) - MV (t H )HB (t)] - p H B - HB (t)), bB [ξ(m)ρBMV (t p - τB )HB (t - τ B )B(t - τB ) - MV (t)HB (t)B(t )] + αB (B0 bP ξ(m)ρP MV (t p - τP )H B (t - τP )B(t - τP )+ αP (P 0 - P (t)), § d § dt HB (t) = § bHB MV (t)HB (t)B(t)+ αB (H0 d B(t) = dt d P (t) = dt d dt F (t) = ρF P (t) - γFV F (t)Vf (t) - αF F (t), - B(t)), где ξ(m) = 1 - m/C0. Биологический смысл параметров системы пояснен в табл. 1 и 2. Для определения решения системы при t> 0 достаточно задать значения MV (t) при -τ5 :( t :( 0, где τ5 = max{τ E,τ B, τE, τB, τP }, значения HE (t) при - max{τ E, τE } :( t :( 0, значения E(t) при -τE :( H H H H t :( 0, значения HB (t) при - max{τ B , τB , τP } :( t :( 0, значения B(t) при - max{τB , τP } :( t :( 0 и ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ РЕШЕНИЙ СИСТЕМ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ 689 Параметр Биологический смысл параметра αM αE H αB H αE αB αP αF τ E H τ B H τE τB τP ρE H ρB H ρE ρB ρP bB p bP p bE H bB H bE p bHE p bHB p ρF bCE bEC Константа скорости потери макрофагом стимулированного состояния Константа скорости потери стимулированного состояния T-хелперами Th1 Константа скорости потери стимулированного состояния T-хелперами Th2 Константа скорости естественной гибели цитотоксических T-лимфоцитов-эффекторов Константа скорости естественной гибели B-лимфоцитов Константа скорости естественной гибели плазматических клеток Константа скорости естественной гибели антител Продолжительность цикла деления Т-хелперов Th1 Продолжительность цикла деления Т-хелперов Th2 Продолжительность цикла делений Т-лимфоцитов-эффекторов Продолжительность цикла делений B-лимфоцитов Продолжительность цикла делений и дифференцировки B-лимфоцитов до появления плазматических клеток Число потомков T-хелпера Th1 в результате 1-го цикла деления Число потомков T-хелпера Th2 в результате 1-го цикла деления Число потомков Т-лимфоцита-эффектора в результате 1-го цикла делений Число потомков B-лимфоцита в результате 1-го цикла делений Число плазматических клеток-потомков В-лимфоцита в результате одного цикла делений Константа скорости стимуляции B-лимфоцита при описании числа B-лимфоцитов Константа скорости стимуляции B-лимфоцита при описании числа плазматических клеток Константа скорости стимуляции Т-хелпера Th1 Константа скорости стимуляции Т-хелпера Th2 Константа скорости стимуляции Т-лимфоцита-эффектора Коэффициент, характеризующий расход Т-хелперов Th1 на стимуляцию Т-лимфоцитов-эффекторов Коэффициент, характеризующий расход Т-хелперов Th2 на стимуляцию B-лимфоцитов Константа скорости синтеза молекул IgG одной плазматической клеткой Константа скорости разрушения гепатоцитов Т-лимфоцитами-эффекторами Константа скорости гибели Т-лимфоцитов-эффекторов вследствие разрушения зараженных клеток Таблица 1. Параметры развития иммунного ответа. значения остальных переменных при t = 0. Однако для единообразия далее мы будем предполагать, что начальные значения всех переменных заданы при -τ5 :( t :( 0. Обозначив вектор переменных системы (2.1) через u(t)= (Vf (t), CV (t), m(t), MV (t), HE (t), E(t), HB (t), B(t),P (t),F (t))T , а вектор параметров через p = (M 0,H0 ,H0 , E0, B0,P 0,C0, αM , αE , αB , αE , αB , αP , αF , ρE , ρB , ρE , ρB , ρP , ρF , bE , bB , E B bE B P H H H H H H HE HB T p , bp , bp , γMV , γFV , σ, bCE , bm, αm, ν, n, γVC , γVM , γVF , bp , bp , bEC ) , 690 М. Ю. ХРИСТИЧЕНКО и др. Parameter The biological meaning of the parameter αM αE H αB H αE αB αP αF τ E H τ B H τE τB τP ρE H ρB H ρE ρB ρP bB p bP p bE H bB H bE p bHE p bHB p ρF bCE bEC Rate constant for the loss of the stimulated state by the macrophage Rate constant of loss of stimulated state by Th1 T-helpers Rate constant of loss of stimulated state by Th2 T-helpers Rate constant of natural death of cytotoxic T-lymphocytes-e ectors Rate constant of natural death of B-lymphocytes Rate constant of natural death of plasma cells Rate constant of natural death of antibodies Duration of the division cycle of Th1 T-helpers Duration of the division cycle of Th2 T-helpers Duration of the division cycle of T-lymphocyte-e ectors Duration of the division cycle of B-lymphocytes Duration of the cycle of divisions and di erentiation of B-lymphocytes before the appearance of plasma cells Number of children of a Th1 T-helper as a result of the 1th division cycle Number of children of a Th2 T-helper as a result of the 1th division cycle Number of children of a T-lymphocyte-e ector as a result of the 1th division cycle Number of children of a B-lymphocyte as a result of the 1th division cycle Number of plasma cells-descendants of a B-lymphocyte as a result of one cycle of divisions B-lymphocyte stimulation rate constant in describing the number of B-lymphocytes B-lymphocyte stimulation rate constant in describing the number of plasma cells Th1 T-helper stimulation rate constant Th1 T-helper stimulation rate constant T-lymphocyte-e ector stimulation rate constant Coe cient characterizing the consumption of Th1 T-helpers for stimulation of T-lymphocytes-e ectors Coe cient characterizing the consumption of Th2 T-helpers for stimulation of T-lymphocytes-e ectors Rate constant of synthesis of IgG molecules by one plasma cell Rate constant of destruction of hepatocytes by T-lymphocytes-e ectors Rate constant of death of T-lymphocyte-e ectors due to the destruction of infected cells Table 1. Parameters of immune response development. рассматриваемую модель можно записать в виде du dt (t)= F (u(t), u(t - τ1),... , u(t - τ5), p) , (2.2) где F : R10×(6) → R10 - рациональная 10-компонентная вектор-функция, а 0 < τ1 < ... < τ5 - задержки. Задача Коши для системы (2.1) с неотрицательными начальными значениями и неотрицательными параметрами глобально разрешима на любом конечном временном интервале. Данное утверждение можно доказать методом шагов Беллмана, рассматривая линейную систему обыкновенных дифференциальных уравнений, мажорирующую правую часть исходной системы, см. [1, 13]. Кроме того, в этих работах было доказано, что если дополнительно выполнены неравенства γVM > 0 и CV (0) + m(0) :( C0, то компоненты Vf (t), MV (t), CV (t), m(t) решения задачи Коши для этой системы ограниченны глобально, причем CV (t)+ m(t) :( C0 при всех t> 0. ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ РЕШЕНИЙ СИСТЕМ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ 691 Параметр Биологический смысл параметра bm αm M 0 H0 E H0 B E0 B0 P 0 C0 γMV γVM γFV γVF σ ν n γVC Константа скорости разрушения зараженных гепатоцитов вследствие цитопатичности вирусов Константа скорости регенерации гепатоцитов Концентрация Ia-несущих макрофагов в лимфоузле Концентрация специфических Т-хелперов Th1 в лимфоузле Концентрация специфических Т-хелперов Th2 в лимфоузле Концентрация специфических предшественников для Т-лимфоцитов-эффекторов в лимфоузле Концентрация специфических B-лимфоцитов в лимфоузле Концентрация специфических плазматических клеток в лимфоузле Концентрация гепатоцитов в печени Константа скорости антигенной стимуляции макрофагов в лимфоузле Константа скорости связывания антигенных частиц макрофагами лимфоузла Константа скорости связывания 1 молекулы IgG с частицей HBsAg Константа скорости нейтрализации вируса гепатита молекулами IgG Константа скорости заражения гепатоцитов Константа скорости секреции частиц HBsAg одним гепатоцитом в сутки Количество частиц HBsAg, высвобождающееся при разрушении гепатоцита Т-лимфоцитом-эффектором Константа скорости адсорбции вирусов незараженными клетками органа-мишени Таблица 2. Параметры иммунного гомеостаза и развития вирусной инфекции в клетках-мишенях. Parameter The biological meaning of the parameter bm αm M 0 H0 E H0 B E0 B0 P 0 C0 γMV γVM γFV γVF σ ν n γVC Rate constant of destruction of infected hepatocytes due to the cytopathicity of viruses Hepatocyte regeneration rate constant Concentration of Ia-carrier macrophages in a lymph node Concentration of speci c Th1 T-helpers in the lymph node Concentration of speci c Th2 T-helpers in the lymph node Concentration of speci c precursors for T-lymphocyte-e ectors in a lymph node Concentration of speci c B-lymphocytes in a lymph node Concentration of speci c plasma cells in a lymph node Concentration of hepatocytes in the liver Rate constant of antigenic stimulation of macrophages in a lymph node Rate constant of binding of antigenic particles by macrophages of a lymph node Binding rate constant of 1 IgG molecule to HBsAg particle Rate constant of neutralization of the hepatitis virus by IgG molecules Rate constant of hepatocyte infection Rate constant of secretion of HBsAg particles by one hepatocyte per day Number of HBsAg particles released during the destruction of a hepatocyte by a T-lymphocyte-e ector Virus adsorption rate constant by uninfected cells of a target organ Table 2. Parameters of immune homeostasis and development of viral infection in target cells. 692 М. Ю. ХРИСТИЧЕНКО и др. 3. Вычисление и анализ стационарных решений Неотрицательные стационарные решения системы вида (2.2) при заданных значениях параметров являются неотрицательными решениями системы алгебраических уравнений G(u, p)= F(u,... ,u, p)= 0. (3.1) 6 Для численного решения систем алгебраических уравнений общего вида в настоящее время используют методы вычислительной алгебры (см. [9]), которые также называют методами символьных вычислений. В частности, такие методы реализованы в процедуре NSolve пакета Mathematica. Эта процедура основана на вычислении базиса Гребнера с использованием мономиального упорядочения и теоретически позволяет решить произвольную систему конечного числа алгебраических уравнений. Однако с ростом суммарной степени рассматриваемых уравнений сложность вычислений растет, вообще говоря, быстрее, чем экспоненциально, и для успешного нахождения решений могут потребоваться предварительные упрощающие преобразования. Для модели (2.1) система (3.1) имеет следующий вид: νCV + nbCECV E - γVF Vf F - γVM Vf - γVCVf fC0 - CV - ml = 0, σVf fC0 - CV - ml - bCECV E - bmCV = 0, bCECV E + bmCV - αmm = 0, γMV M 0Vf - αM MV = 0, [bE (ξ(m)ρE - 1)]MV HE - bHE MV HEE + αE (H0 - HE )= 0, H H p H E (3.2) p [bE (ξ(m)ρE - 1)]MV HEE - bECCV E + αE (E0 - E)= 0, [bB (ξ(m)ρB - 1)]MV HB - bHB MV HBB + αB (H0 - HB )= 0, H H p H B p p [bB (ξ(m)ρB - 1)]MV HBB + αB (B0 - B)= 0, [bP ξ(m)ρP ]MV HBB + αP (P 0 - P )= 0, ρF P - γFV F Vf - αF F = 0, где ξ(m)= 1 - m/C0. В работе [12] показано, что решить систему (3.2) при фиксированных значениях параметров можно с помощью процедуры NSolve пакета Mathematica, выполнив предварительно упрощающие преобразования, предложенные в работе [18] и сводящие эту систему к системе из четырех уравнений, а именно, к следующей системе уравнений относительно переменных Vf , CV , HB и HE : [bB (ξ(m)ρB - 1)]MV HB - bHB R + αB (H0 - HB )= 0, H H p H B MV HBB - R = 0, [bE (ξ(m)ρE - 1)]MV HE - bHE MV HEE + αE (H0 - HE )= 0, H H p H E где p [bE (ξ(m)ρE - 1)]MV HEE - bECCV E + αE (E0 - E)= 0, m = σVf (C0 - CV ) αm + σVf , E = αmm - bmCV bCECV , MV = γMV M 0 Vf , αM ( γFV Vf + αF \ αP R = (νCV + nbCECV E - γVM Vf - γVCVf (C0 - CV - m)) ρ γ V - P0 bP , F bP B VF f p ξ(m)ρP p B - P = p ξ(m)ρP R + P 0, F = ρF P [b (ξ(m)ρ 1)]R , B = + B0 αP γFV Vf + αF αB - функции, полученные посредством явного выражения остальных переменных системы (3.2) через переменные Vf , CV . ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ РЕШЕНИЙ СИСТЕМ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ 693 Для трассирования решений системы (3.1) по параметрам можно использовать описанный ниже алгоритм, являющийся обобщением алгоритма, предложенного в [16] для трассирования стационарных решений модели динамики инфекции, вызванной вирусами лимфоцитарного хориоменингита (ВЛХМ). Пусть в пространстве значений вектора параметров p задана некоторая гладкая рационально параметризуемая кривая, не имеющая самопересечений. Пусть требуется найти решения системы (3.1) как функции параметра этой кривой. Будем обозначать параметр через s, а интервал его варьирования через η, предполагая этот интервал замкнутым. Введем обозначение S(s, u)= G(u, p(s)). Поскольку S(s, u) является рациональной функцией переменных s и u, множество A = {(s, u): S(s, u)= 0, s ∈ η, u 0} является частью некоторой действительной алгебраической кривой, см. [11]. Поставленная задача сводится к разбиению множества A на подмножества вида C = {(s, u(s)) : s ∈ ηC }, (3.3) где u(s) - однозначная аналитическая функция переменного s, определенная в некотором интервале ηC ⊂ η, и точки, в которых выполняется равенство det ( ∂S ∂u \ (s, u) = 0. (3.4) Отметим, что поскольку A является частью алгебраической кривой, число таких точек конечно. u u u i i+1 1 1 u u i i+1 2 2 i u i+1 u3 3 3 u~ i+1 si si+1 s Рис. 1. i-й шаг второго этапа трассирования решений Fig. 1. i-th step of the second stage of solution tracing Алгоритм построения подмножеств (3.3) состоит из двух этапов. На первом этапе задается достаточно мелкая сетка в интервале η с узлами s1 < ... < sN , где s1 - начало, а sN - конец этого интервала, и вычисляются множества решений системы (3.1) вида j U(si)= {ui : j = 1,... , mi} (3.5) с помощью процедуры NSolve пакета Mathematica, где mi - число найденных неотрицательных решений в узле si. На втором этапе выполняется соединение элементов множеств (3.5), полученных для соседних узлов. Этот этап состоит из N -1 шага. На i-ом шаге рассматривается интервал [si, si+1]. Сначала стационарные решения из множества U(si) трассируются по параметру от значения si до значения si+1 путем численного решения задач Коши u(si)= ui , du = - ( ∂S \-1 ∂S , j = 1,... ,m , (3.6) j ds ∂u ∂s i с помощью стандартной процедуры численного интегрирования ode45 среды MATLAB, основанной на методах Рунге-Кутта четвертого и пятого порядков. Интегрирование останавливается, 694 М. Ю. ХРИСТИЧЕНКО и др. если переменная s достигает заданного финального значения si+1, либо если в процессе интегрирования обнаружится, что функция f (s) = det(∂S/∂u(s, u)) обращается в нуль на траектории интегрирования. Отметим, что сетку интегрирования по s процедура ode45 выбирает самостоятельно. В результате трассирования по параметру в узле si+1 мы получаем множество U (si+1)= {ui+1 : j = 1,... , mi+1}, j mi+1 где :( m i+1 . После этого элементы множества U (si+1) уточняются элементами множества U(si+1) путем выбора элементов множества U(si+1), наиболее близких по норме к элементам множества mi+1 U (si+1). Если < mi+1 , то оставшиеся m i+1 - mi+1 стационарных решения, не полученные в результате трассирования по параметру от значения si до значения s i+1 , необходимо трассировать по параметру от значения si+1 до значения si путем решения задач Коши, аналогичных задачам Коши (3.6). Рисунок 1 поясняет i-й шаг второго этапа. Круглыми точками обозначены элементы множеств U(si) и U(si+1), квадратными - конечные точки численного интегрирования. Черные линии обозначают трассирование по параметру s путем численного интегрирования и уточнение, красные - линейные аппроксимации искомых подмножеств вида (3.3) в интервале [si, si+1]. Анализ устойчивости найденного стационарного состояния системы с запаздыванием при фиксированных значениях параметров сводится к решению нелинейной проблемы собственных значений с полиномиальным и экспоненциальным вхождением спектрального параметра. Эффективный метод решения таких проблем описан в работах [16, 18]. 4. Результаты численных экспериментов С помощью предложенной технологии была исследована зависимость стационарных решений модели Марчука-Петрова от константы скорости антигенной стимуляции макрофагов в лимфоузле γMV . Этот параметр варьировался в интервале [1,7 · 10-13, 1,7 · 10-11]. Вычисление стационарных решений и анализ их зависимости от параметра γMV выполнялся при трех различных наборах параметров, которые мы в дальнейшем будем называть наборами (a), (b) и (c). Набор p (a) соответствует острому течению гепатита B. Значения параметров из этого набора приведены в табл. 3. Набор (b) состоит из тех же значений параметров, что и набор (a), за исключением параметра ν, который был уменьшен в 3 раза. Набор (c) состоит из тех же значений параметров, что и набор (b), за исключением параметров ρF , bm, bCE , bE , которые были выбраны равными 8,5·106, 5,2·10-2, 5,2·10-5, 4,1·10-8 соответственно. Биологически, набор значений параметров (b) соответствует вирусной инфекции с малой скоростью размножения вирусов, что является одной из причин хронизации инфекции [13]. Набор значений параметров (c) отвечает ситуации хронической инфекции с более эффективной стимуляцией цитотоксического T-клеточного иммунного ответа. В дополнение к варьированию значений параметров модели, мы варьировали чувствительность иммунного ответа к поражению органа мишени, которая описывается с помощью функции ξ(m). Так, в экспериментах с набором (c) вместо функции ξ(m) = max{1 - m/C0, 0} использовались функции ξ(m) = max{(1 - m/C0)5, 0} и ξ(m) = max{(1 - m/(0,2C0)), 0}. На рис. 3-4 приведены зависимости стационарных решений модели Марчука-Петрова от параметра γMV для наборов (a) и (b). На рис. 5 приведена зависимость стационарных решений модели Марчука-Петрова от параметра γMV в интервале [1,5 · 10-11, 1,7 · 10-11] для набора (с) при ξ(m)= max{(1 - m/C0)5, 0}. На рис. 6 приведена зависимость стационарных решений модели Марчука-Петрова от параметра γMV в интервале [1,6 · 10-11, 2,5 · 10-11] для набора (с) при ξ(m)= max{(1 - m/(0,2C0)), 0}. Выбор альтернативных параметризаций функции ξ(m) обусловлен следующими соображениями. Поражение клеток печени вследствие цитопатических эффектов вируса, а также вследствие их уничтожения цитотоксическими лимфоцитами, вызывает супрессию антиген-специфического иммунного ответа. Данная обратная связь параметризована с помощью функции ξ(m), где m - количество пораженных клеток печени, m ∈ [0,C0]. В исходной модели параметризация имеет вид ξ(m) = 1 - m/C0. Функция ξ(m) принимает значения от 0 (полная супрессия иммунного ответа при полном разрушении печени) до 1 (полноценный иммунный ответ при отсутствии повреждений печени). В данной работе предложены две альтернативные параметризации этой ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ РЕШЕНИЙ СИСТЕМ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ 695 Параметр Оценка Параметр Оценка αM αE H αB H αE αB αP αF τ E H τ B H τE τB τP ρE B H , ρH ρE, ρB ρP B bp P bp bE H bB H E bp HE bp 1,2 сут-1 1,0 сут-1 1,0 сут-1 0,4 сут-1 0,1 сут-1 0,4 сут-1 0,043 сут-1 0,6 сут 0,6 сут 2,0 сут 2,0 сут 3,0 сут 2 16 3 2,2 · 10-9 (мл/клеток)2 · сут-1 4,7 · 10-12 (мл/клеток)2 · сут-1 4,5 · 10-5 (мл/клеток) · сут-1 4,5 · 10-5 (мл/клеток) · сут-1 1,5 · 10-8 (мл/клеток)2 · сут-1 1,5 · 10-14 (мл/клеток)2 · сут-1 bHB p ρF bCE bEC bm αm M 0 H0 E H0 B E0 B0 P 0 C0 γMV γVM γFV γVF σ ν n γVC 2,2 · 10-13 (мл/клеток)2 · сут-1 1,7 · 108(частиц/клеток) · сут-1 1,1 · 10-6 (мл/частиц) · сут-1 2,7 · 10-7 (мл/клеток) · сут-1 0,068 сут-1 0,15 сут-1 6,0 · 105 клеток/мл 6,0 · 102 клеток/мл 6,0 · 101 клеток/мл 6,0 · 102 клеток/мл 6,0 · 102 клеток/мл 2,6 · 10-1 клеток/мл 3,0 · 108 клеток/мл 1,6 · 10-11 (мл/клеток) · сут-1 0,4 сут-1 1,4 · 10-9 (мл/частиц) · сут-1 5,0 · 10-10 (мл/частиц) · сут-1 3,8 · 10-12 (мл/клеток) · сут-1 83 сут-1 5 4,2 · 10-14 (мл/клеток) · сут-1 Таблица 3. Значения параметров модели (2.1), соответствующие острому гепатиту B (см. [15]). Parameter Estimate Parameter Estimate αM αE H αB H αE αB αP αF τ E H τ B H τE τB τP ρE B H, ρH ρE, ρB ρP B bp P bp bE H bB H E bp HE bp 1.2 day-1 1.0 day-1 1.0 day-1 0.4 day-1 0.1 day-1 0.4 day-1 0.043 day-1 0.6 day 0.6 day 2.0 day 2.0 day 3.0 day 2 16 3 2.2 · 10-9 (ml/cells)2 · day-1 4.7 · 10-12 (ml/cells)2 · day-1 4.5 · 10-5 (ml/cells) · day-1 4.5 · 10-5 (ml/cells) · day-1 1.5 · 10-8 (ml/cells)2 · day-1 1.5 · 10-14 (ml/cells)2 · day-1 bHB p ρF bCE bEC bm αm M 0 H0 E H0 B E0 B0 P 0 C0 γMV γVM γFV γVF σ ν n γVC 2.2 · 10-13 (ml/cells)2 · day-1 1.7 · 108(particles/cells) · day-1 1.1 · 10-6 (ml/particles) · day-1 2.7 · 10-7 (ml/cells) · day-1 0.068 day-1 0.15 day-1 6.0 · 105 cells/ml 6.0 · 102 cells/ml 6.0 · 101 cells/ml 6.0 · 102 cells/ml 6.0 · 102 cells/ml 2.6 · 10-1 cells/ml 3.0 · 108 cells/ml 1.6 · 10-11 (ml/cells) · day-1 0.4 day-1 1.4 · 10-9 (ml/particles) · day-1 5.0 · 10-10 (ml/particles) · day-1 3.8 · 10-12 (ml/cells) · day-1 83 day-1 5 4.2 · 10-14 (ml/cells) · day-1 Table 3. Model parameter values (2.1) corresponding to acute hepatitis B (see [15]). 696 М. Ю. ХРИСТИЧЕНКО и др. функции: ( m \k а также ξ(m) ≡ ξk(m)= ⎧ 1 - C0 m , (4.1) ⎨1 - ξ(m) ≡ ξr (m)= rC0 , m :( rC0, (4.2) ⎩0, m > rC0, которые эквивалентны исходной функции при параметрах k =1 и r = 1. При увеличении k > 1 или уменьшении r < 1 функция ξ(m) быстрее убывает от единицы до нуля (см. рис. 2), что соответствует более сильной чувствительности супрессии иммунного ответа к степени поражения печени. Альтернативные параметризации функции ξ(m) позволили выявить различные бистабильные режимы со свойством гистерезиса (рис. 5 и 6). Рис. 2. Параметризация эффекта супрессии антиген-специфического иммунного ответа при повреждении печени с помощью функций ξk (m) и ξr (m) (уравнения (4.1), (4.2)). Данные зависимости показаны красным цветом слева и справа, соответственно. Синяя линия соответствует зависимости ξ(m) в исходной формулировке модели Марчука-Петрова [1]. Fig. 2. Parameterization of the antigen-speci c immune response suppression e ect in liver damage using the ξk (m) and ξr(m) functions (Eqs. (4.1) and (4.2)). These dependencies are shown in red on the left and right, respectively. The blue line corresponds to the dependence of ξ(m) in the original formulation of the Marchuk- Petrov model [1]. 1. Моностабильность. Трассировка стационарных решений по параметру γMV , представленная на рис. 3, показывает, что для возникновения условия хронической инфекции, т. е. устойчивого сосуществования популяции вирусов и иммунного ответа, требуется существенное снижение эффективности антигенпрезентирования до уровня γMV ≈ 3 · 10-12. При этом уменьшение ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ РЕШЕНИЙ СИСТЕМ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ 697 скорости размножения вирусов ν в 3 раза существенно расширяет интервал значений параметра γMV , в котором возможно формирование стационарных решений, соответствующих хроническому течению вирусного гепатита В с высокой вирусной нагрузкой (рис. 4). Стоит отметить, что результат, проиллюстрированный на рис. 4, соответствует полученным ранее теоретическим результатам. В работе [4] было приведено следующее условие асимптотической устойчивости тривиального стационарного решения с Vf = 0: ( G = γVF ρF P 0 αF + γVM + γVCC 0\ bCEE0 + bm - σC 0 ν + nbCEE0 > 0. Нетрудно проверить, что это условие выполняется для набора параметров (b) и не выполняется для набора (a), что соответствует результатам, проиллюстрированным на рис. 3-4. 2. Бистабильность. Бистабильность является свойством динамической системы, предполагающим сосуществование двух устойчивых положений равновесия в фазовом пространстве системы. Для модели Марчука-Петрова они соответствуют тяжелой и легкой формам заболевания, если рассматривать вирусную нагрузку как индикатор тяжести вирусного гепатита В. Наличие бистабильности возможно в определённых областях параметров модели, поиск которых был осуществлен с использованием разработанных нами методов. На рис. 5, соответствующем набору параметров (с), показано сосуществование трех нетривиальных положений равновесия, два из которых являются одновременно устойчивыми для некоторого интервала значений параметра γMV - константы скорости активации антигенпрезентирующих клеток. Ширина интервала сосуществования двух устойчивых стационарных решений достаточно мала, что определяет сложность поиска таких режимов. Чувствительность иммунного ответа к повреждению клеток печени, описываемая с помощью функции ξ(m), существенно влияет на размер области бистабильности по параметру γMV . Увеличение чувствительности может приводить к появлению свойства ультрачувствительности, которое важно для реализации бистабильности. На рис. 6 показано, что переход от ξ(m) = max{(1 - m/C0)5, 0} к ξ(m) = max{(1 - m/(0,2C0)), 0}, приводящий к увеличению производной данной функции, приводит к сужению области существования устойчивого решения с меньшей вирусной нагрузкой. Отметим, что бистабильность является частным случаем свойства мультистабильности динамических систем. Известно, что реализация свойства мультистабильности является нетривиальной задачей и определяется регуляторными блоками математических моделей, см. [19]. Результаты анализа множества стационарных точек в модели Марчука-Петрова, представленные на рис. 5-6, демонстрируют наличие 4-х стационарных решений, из которых два являются неустойчивыми. Поиск областей значений параметров, в которых возможно сосуществование более чем двух устойчивых стационарных решений, требует дальнейших систематических исследований в данной многопараметрической модели. 3. Гистерезис. Характер зависимости бистабильности от параметра γMV , продемонстрированный на рис. 5-6, позволил выявить свойство гистерезиса в модели Марчука-Петрова. Гистерезис является свойством динамической системы, которое определяет зависимость её состояния при изменении параметров от предыстории. Так, если система находилась в состоянии с более высокой вирусной нагрузкой, как показано на рис. 5 (верхняя ветвь), то увеличение параметра скорости активации макрофагов от значения 1,5 · 10-11 до 1,52 · 10-11 не приведет к переходу в состояние с более низкой вирусной нагрузкой. Для этого потребуется его увеличение до значения 1,58 · 10-11 , при котором происходит потеря устойчивости и переход системы в другое стационарное состояние с более низкой вирусной нагрузкой. Соответственно, если имеет место снижение скорости активации макрофагов до значения 1,52 · 10-11, то система остается на нижней ветви устойчивого стационарного состояния. Случайные флуктуации данного параметра в области бистабильности будут проявляться в неоднозначности траектории системы, что связано с тем, на какой из ветвей находилась система в процессе изменений. 4. Бифуркации. Зависимости стационарных решений от параметра γMV , изображенные на рис. 3, 5 и 6, имеют точки бифуркации. Они хорошо видны на графиках зависимости переменной E от этого параметра. На рис. 3 видны две точки бифуркации. С ростом значения параметра при 698 М. Ю. ХРИСТИЧЕНКО и др. γMV ≈ 3 · 10-12 происходит потеря устойчивости стационарного решения, изображенного красной линией, при переходе комплексно-сопряженной пары ведущих (с максимальной вещественной частью) собственных значений линеаризованной относительно этого стационарного решения системы из левой полуплоскости в правую. Наши дополнительные расчеты, результаты которых не вошли в данную работу, показали, что эта бифуркация приводит к образованию в окрестности неустойчивого стационарного решения устойчивого периодическое решения, период и размах которого растут с ростом значения параметра. Таким образом, в этой точке происходит бифуркация Андронова-Хопфа, см. [17]. Вторая бифуркации происходит при уменьшении значении параметра в точке γMV ≈ 8,2 · 10-12. Неустойчивые стационарные решения, изображенные зеленой и синей линиями, сходятся в одну точку и на этом обрываются. Такая точка бифуркации называется точкой поворота, см. [17]. Ведущие собственные значения линеаризованной относительно стационарного решения в этой точке системы представляют собой комплексно-сопряженную пару с положительной вещественной частью. На рис. 5 видны две точки бифуркации: при γMV ≈ 1,52 · 10-11 соединяются устойчивое стационарное решение, изображенное синей линией, и неустойчивое стационарное решение, изображенное зеленой линией, а при γMV ≈ 1,59 · 10-11 соединяются устойчивое стационарное решение, изображенное красной линией, и неустойчивое стационарное решение, изображенное зеленой линией. Это так называемые точки касательной бифуркации, см. [17]. В этих точках ведущее собственное значение линеаризованной относительно стационарного решения системы вещественно и равно нулю. На рис. 6 видны три точки бифуркации: при γMV ≈ 1,62·10-11 сходятся устойчивое стационарное решение, изображенное синей линией, и неустойчивое стационарное решение, изображенное зеленой линией, а при γMV ≈ 2,25 · 10-11 сходятся устойчивое стационарное решение, изображенное красной линией, и неустойчивое стационарное решение, изображенное зеленой линией. Эти точки являются точками касательной бифуркации. При γMV ≈ 1,7 · 10-11 происходит потеря устойчивости стационарного решения, изображенного красной линией. Так же как в аналогичном случае, изображенном на рис. 3, эта потеря устойчивости происходит при переходе комплексносопряженной пары ведущих собственных значений из левой полуплоскости в правую и приводит к образованию в окрестности неустойчивого стационарного решения устойчивого периодического решения, период и размах которого растут с ростом значения параметра. Таким образом, в этой точке происходит бифуркация Андронова-Хопфа. 5. Выводы В данной работе предложены эффективные методы расчета и численного анализа устойчивости стационарных решений нелинейных математических моделей на основе дифференциальных уравнений с запаздывающим аргументом. С помощью предложенных методов для математической модели противовирусного иммунного ответа Марчука-Петрова впервые проведено исследование стационарных решений, соответствующих хроническим формам течения инфекции вирусами гепатита В. Было показано наличие бистабильности для данной модели и свойство гистерезиса, а также определены области в пространстве параметров модели, в которых данные свойства имеют место. Наличие бистабильности позволяет предлагать различные подходы к лечению неблагоприятных вариантов хронического гепатита В, в частности, переводить систему в состояние с более низкой вирусной нагрузкой, на основе разработанного нами ранее метода оптимальных возмущений, см. [5, 6]. Вместе с тем, хроническое течение вирусного гепатита В может иметь место в виде осциллирующей динамики. Исследование соответствующих решений и разработка эффективных методов бифуркационного анализа для систем функционально-дифференциальных уравнений определяет одно из актуальных направлений дальнейших исследований, с приложениями в области патогенеза социально значимых заболеваний, таких как COVID-19, ВИЧ и гепатит С. Финансирование. Разработка алгоритмов и численные эксперименты выполнены при финансовой поддержке Российского научного фонда (проект № 22-11-00025), биологическая интерпретация результатов выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 20-01-00352). ЧИСЛЕННЫЙ АНАЛИЗ СТАЦИОНАРНЫХ РЕШЕНИЙ СИСТЕМ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ 699 1011 109 108 105 2000 1010 109 108 108 107 107 106 104 1800 1600 107 106 105 106 105 104 105 104 3 103 102 1400 1200 104 103 102 103 102 10 102 1 101 1000 800 1 101 -10 5 10 15 10-12 101 1 0 -1 10 1 0 -1 5 10 15 10-12 1 0 5 10 15 10-12 5 10 15 10-12 600 5 10 15 10-12 105 104 103 180 160 140 105 60 109 50 40 108 30 102 101 1 0 5 10 15 10-12 120 100 80 60 5 10 15 10-12 104 103 20 10 0 5 10 15 10-12 5 10 15 10-12 107 106 5 10 15 10-12 Рис. 3. Зависимость стационарных решений от параметра γMV ∈ [1,7 · 10-13, 1,7 · 10-11] (сплошная линия - устойчивое решение, пунктирная - неустойчивое) для набора параметров (a), соответствующего острой форме вирусного гепатита В. Fig. 3. Dependence of stationary solutions on the parameter γMV ∈ [1.7 · 10-13, 1.7 · 10-11] (solid line: stable solution, dotted line: unstable) for the set of parameters (a) corresponding to the acute form of viral hepatitis B. 109 108 107 106 105 104 103 102 107 106 105 104 103 102 1 107 106 105 104 103 102 101 104 103 102 101 1400 1200 1000 800 101 1 0 10 1 0 5 10 15 10-12 1 0 -1 5 10 15 10-12 1 0 5 10 15 10-12 5 10 15 10-12 600 5 10 15 10-12 700 600 500 400 300 200 100 0 5 10 15 10-12 140 120 100 80 60 5 10 15 10-12 1400 1200 1000 800 600 5 10 15 10-12 0.36 0.34 0.32 0.3 0.28 0.26 0.24 5 10 15 10-12 109 108 107 5 10 15 10-12 Fig. 4. Зависимость стационарных решений от параметра γMV ∈ [1,7 · 10-13, 1,7 · 10-11] (сплошная линия - устойчивое решение, пунктирная - неустойчивое) для набора параметров (b), в котором скорость размножения вирусов ν уменьшена в 3 раза по сравнению с острой формой, что расширяет интервал значений параметра γMV , в котором возможно формирование стационарных решений, соответствующих хроническому течению вирусного гепатита В с высокой вирусной нагрузкой. Fig. 4. Dependence of stationary solutions on the parameter γMV ∈ [1.7 · 10-13, 1.7 · 10-11] (solid line: stable solution, dotted line: unstable) for the set of parameters (b), in which the rate of reproduction of viruses ν is reduced by a factor of 3 compared to the acute form, which expands the range of values of the parameter γMV in which the formation of stationary solutions is possible corresponding to the chronic course of viral hepatitis B with a high viral load. 700 М. Ю. ХРИСТИЧЕНКО и др. 5252 52. 529 528 529 528 520 7 5422 529 528 52- 52- 520 52- 520 52 526 5522 5222 520 527 526 527 526 527 526 524 . 22 922 524 5 525 352 510 5100 5151-0 52355 524 525 5 2 35 510 5100 5151-0 52355 524 525 5 2 35 510 5100 5151-0 52355 525 5 2 35 510 5100 5151-0 52355 822 22 510 5100 5151-0 52355 527 526 524 525 510 5100 5151-0 52355 562 542 552 522 . 2 92 82 -2 510 5100 5151-0 52355 527 526 510 5100 5151-0 52355 51- 517 514 5 219 21- 217 214 510 5100 5151-0 52355 529 528 52- 520 510 5100 5151-0 52355 Рис. 5. Зависимость стационарных решений от параметра γMV ∈ [1,5 · 10-11, 1,7 · 10-11] (сплошная линия - устойчивое решение, пунктирная - неустойчивое) для набора параметров (c), который отвечает хронической инфекции с более эффективной стимуляцией цитотоксического T-клеточного иммунного ответа, при ξ(m) = (1 - m/C0)5. Демонстрирует появление бистабильности и гистерезиса в модели Марчука-Петрова. Fig. 5. Dependence of stationary solutions on the parameter γMV ∈ [1.5 · 10-11, 1.7 · 10-11] (solid line: stable solution, dotted line: unstable) for the set of parameters (c), which corresponds to chronic infection with more e ective stimulation of the cytotoxic T-cell immune response, at ξ(m) = (1 - m/C0)5. It demonstrates the bistability and hysteresis in the Marchuk-Petrov model. 54. 8 9 9 54 5544 549 540 540 543 7 5444 . 44 548 542 5 654 510 512 212 213 54655 543 54 542 5 5 4 65 510 512 212 213 54655 543 54 542 5 5 4 65 510 512 212 213 54655 54 54 45 944 5 4 044 5 510 512 212 213 54655 44 510 512 212 213 54655 543 547 554 544 . 4 543 2 518 54- 542 545 510 512 212 213 54655 -4 94 04 510 512 212 213 54655 547 510 512 212 213 54655 5 418 4 510 512 212 213 54655 549 540 510 512 212 213 54655 Fig. 6. Зависимость стационарных решений от параметра γMV ∈ [1,6 · 10-11, 2,5 · 10-11] (сплошная линия - устойчивое решение, пунктирная - неустойчивое) для набора (c), который отвечает хронической инфекции с более эффективной стимуляцией цитотоксического T-клеточного иммунного ответа, при ξ(m)= (1 - m/(0,2C0)). Демонстрирует появление бистабильности и гистерезиса в модели Марчука-Петрова. Fig. 6. Dependence of stationary solutions on the parameter γMV ∈ [1.6 · 10-11, 2.5 · 10-11] (solid line: stable solution, dotted line: unstable) for the set (c), which responds to chronic infection with more e ective stimulation of the cytotoxic T-cell immune response, at ξ(m) = (1 - m/(0.2C0)). It demonstrates the bistability and hysteresis in the Marchuk-Petrov model.

Об авторах

М. Ю. Христиченко

Институт вычислительной математики им. Г. И. Марчука РАН

Автор, ответственный за переписку.
Email: misha.hrist@gmail.com
Москва, Россия

Ю. М. Нечепуренко

Институт вычислительной математики им. Г. И. Марчука РАН

Email: yumnech@yandex.ru
Москва, Россия

Д. С. Гребенников

Институт вычислительной математики им. Г. И. Марчука РАН

Email: dmitry.ew@gmail.com
Москва, Россия

Г. А. Бочаров

Институт вычислительной математики им. Г. И. Марчука РАН

Email: gbocharov@gmail.com
Москва, Россия

Список литературы