CC BY
45
Челябинский физико-математический журнал. 2016. Т. 1, вып. 1. С. 73-80.
УДК 517.929
МОДЕЛИРОВАНИЕ ДИНАМИЧЕСКОГО ПРОЦЕССА С ЭФФЕКТОМ ЗАПАЗДЫВАНИЯ
С. А. Никитина", А. С. Скорынин6
Челябинский государственный университет, Челябинск, Россия "nikitina@csu.ru; 6anton.skorynin@gmail.com
В работе рассматриваются вопросы, связанные с моделированием систем с последействием. При математическом описании динамического процесса обычно предполагают, что его поведение зависит только от текущего состояния. Однако в некоторых случаях поведение системы также зависит от информации о её предыдущих состояниях. При этом появление запаздывания в системе может существенно изменить её динамику. Для моделирования таких систем применяются специальные методы.
Ключевые слова: математическое моделирование, динамическая система с последействием, эффект запаздывания.
Постановка задачи
Опишем модель инфекционного заболевания [1]. Под термином «инфекционная болезнь» понимается отображение взаимоотношений, установившихся между членами биоценоза, один из которых благодаря механизмам патогенности способен существовать в другом, а другой благодаря механизмам защиты — оказать противодействие патогенному действию. Одним из основных способов защиты организма от инфекции является система иммунитета.
Суть иммунной реакции на вторжение генетически чужеродного материала (антигена), в том числе и возбудителя болезни, состоит в производстве специфических объектов (антител, клеток-киллеров), способных нейтрализовать или уничтожать антигены. Исходя из этого, инфекционное заболевание рассматривается как конфликт между популяцией возбудителей болезни и иммунной системой организма.
В рамках модели процесс заболевания описывается следующим образом. В момент ¿0 в организм проникает начальная популяция вирусов У0. Вирусы начинают размножаться в клетках органа-мишени и тем самым поражают его. Часть вирусов поступает в кровь, где взаимодействует с рецепторами иммунокомпетентных клеток, в результате происходит стимуляция иммунной системы. Спустя время т после стимуляции в организме появляются клоны плазматических клеток, вырабатывающие антитела, специфичные к вирусам. Антитела связывают вирусы, и от относительной кинетики этих процессов зависит исход болезни. Если вирусы успевают значительно поразить ткань органа, то ухудшается общее состояние организма и, как следствие, работы иммунной системы. Производство антител падает, и снижается вероятность благоприятного исхода.
Основными факторами инфекционного заболевания являются следующие величины. V(¿) — концентрация вирусов (размножающиеся патогенные антигены); ^(¿) — концентрация антител (субстраты иммунной системы, нейтрализующие вирусы: иммуноглобулины, рецепторы клеток); С(¿) — концентрация плазматических
клеток — популяция носителей и продуцентов антител (иммунокомпетентные клетки и иммуноглобулинопродуценты); m(t) — относительная характеристика поражённого органа (масса, число поражённых клеток). Модель сформулирована в виде системы нелинейных дифференциальных уравнений [1; 2]
i f = (в - im
§ = aF(t - т)V(t - т) - цc(C - C*),
(1)
f = pC - ц + niV)F,
l t = °V - Ц^т.
Первое уравнение описывает изменение числа антигенов в организме. Оно может быть переписано в виде
dV = eVdt - YFVdt.
Первый член в правой части уравнения описывает прирост антигенов dV за интервал времени dt за счёт размножения. Естественно, что он пропорционален V и некоторому числу в, которое называется коэффициентом размножения антигенов. Член -YfFVdt описывает число антигенов, нейтрализуемых антителами F за интервал времени dt. Таким образом, число таких вирусов будет пропорционально как количеству антител в организме, так и числу антигенов; y — коэффициент, связанный с вероятностью нейтрализации антигена антителами при встрече с ним.
Второе уравнение в системе
dC = aF (t - т )V (t - т )dt - цс(С - C*)dt
описывает рост плазматических клеток. Первый член в правой части описывает генерацию плазмоклеток, т — время, в течение которого осуществляется формирование каскада плазматических клеток (время запаздывания), a — коэффициент, учитывающий вероятность встречи антиген — антитело, возбуждение каскадной реакции и число образующихся новых клеток. Второй член описывает уменьшение числа плазматических клеток за счёт старения, цс — коэффициент, равный обратной величине их времени жизни.
При получении третьего уравнения в системе подсчитан баланс числа антител, реагирующих с антигеном:
dF = pCdt - njVFdt - цf Fdt.
Первый член справа pCdt описывает генерацию антител плазматическими клетками за интервал времени dt, p — коэффициент производства антител одной плазматической клеткой. Второй член niFVdt описывает уменьшение числа антител в интервале времени dt за счёт связи с антигенами. В самом деле количество выбывающих за интервал времени dt антигенов за счёт нейтрализации их антителами было равно YFVdt. Если на нейтрализацию одного антигена требуется n антител, то мы и приходим к указанному члену формулы. Третий член ц-f Fdt описывает уменьшение популяции антител за счёт старения, где Цf — коэффициент, обратно пропорциональный времени распада антител.
Построенные уравнения не учитывают ослабления жизнедеятельности организма в ходе заболевания, связанного с уменьшением активности органов, обеспечивающих поставку иммунологического материала: лейкоцитов, лимфоцитов, антител
и т. д., необходимых для борьбы с размножающимися антигенами. Производительность таких органов связана с размерами поражения органа-мишени. Рассмотрим уравнение для относительной характеристики поражения органа-мишени. Пусть М — характеристика здорового органа (масса или площадь), а М' — соответствующая характеристика здоровой части поражённого органа. Рассмотрим величину т = 1 — М'/М. Это относительная характеристика поражения органа-мишени. Для непоражённого органа она равна нулю, а для полностью поражённого — единице. Для этой характеристики рассмотрим четвёртое уравнение в системе (1).
Первый член в правой части этого уравнения характеризует степень поражения органа. Предполагается, что за интервал времени увеличение относительной величины поражённого органа пропорционально количеству антигенов, которое описывается членом аV, где а — некоторая константа, своя для каждого заболевания. Уменьшение этой характеристики происходит за счёт восстановительной деятельности организма. Этот член будет зависеть от т с коэффициентом пропорциональности
К системе уравнений (1) присоединим начальные данные при £ = ¿0. Обычно для уравнений с запаздыванием начальные условия задаются на интервале [¿о — т, ¿0]. Однако в биологическом смысле описываемых процессов до момента заражения ¿0 вирусов в организме не было, т. е. V(¿0) = 0. Далее, когда речь будет идти о начальных условиях для уравнений такого типа, задание их в начальной точке означает, что V(¿) = 0 при ¿ <¿0. Имеем
V (¿0 ) = Vо, С (¿0) = С0, ^ (¿0) = Fо, т(£0 ) = т0. 1. Численные методы реализации математической модели
Математическая модель инфекционного заболевания сформулирована в виде системы дифференциальных уравнений с запаздывающим аргументом. Основной подход в построении методов численного решения дифференциальных уравнений с запаздывающим аргументом состоит в объединении какого-либо метода численного интегрирования для обыкновенных дифференциальных уравнений и метода интерполяции. Для аппроксимации запаздывающих переменных, как правило, используют методы интерполяции в виде сплайнов, многочленов Лагранжа и Эрмита и т. д. Причём порядок интерполяции должен согласовываться с порядком аппроксимации разностного уравнения и быть не ниже последнего [3].
Для численного решения данной системы дифференциальных уравнение были реализованы два метода: метод Рунге — Кутта — Фельберга и метод Адамса — Бэшфорта — Мултона. Подробно эти методы представлены в работах [4; 5]. Остановимся лишь на их основных особенностях, связанных с появлением в системе запаздывания.
1.1. Реализация метода Рунге — Кутта — Фельберга
При компьютерном моделировании [6] был применён метод Рунге — Кутта — Фельберга четвёртого и пятого порядка, адаптированный для уравнений с запаздывающим аргументом. Модификация данного алгоритма связана с использованием интерполяции полиномами Эрмита пятого порядка для аппроксимации запаздывающих переменных.
При использовании метода Рунге — Кутта — Фельберга четвёртого и пятого порядка следует учитывать существование точек разрывов для первых семи производных решения. Возникает необходимость организовать вычислительный процесс
таким образом, чтобы точки разрывов производных решения принадлежали множеству узлов сетки интегрирования и не являлись внутренними точками отрезка, на котором определяется интерполяционный многочлен Эрмита.
1.2. Реализация метода Адамса — Бэшфорта — Мултона
Метод Адамса — Бэшфорта — Мултона принадлежит к классу линейных многошаговых методов. Для аппроксимации запаздывающих переменных используются интерполяционные полиномы корректора, формируемые на основе вектора Норд-сика. В процессе решения точки разрывов первых восьми производных решения автоматически совмещаются с узлами сетки интегрирования.
2. Результаты моделирования
Анализ модели заболевания (1) позволяет судить о поведении решения V(¿) — концентрации антигенов при том или ином наборе коэффициентов.
После проведения ряда численных экспериментов, а также в результате консультаций со специалистами в области иммунологии установлено, что модель (1) описывает динамику тех инфекционных заболеваний, которые могут вызывать поражения только одного органа человека. Одним из таких заболеваний является вирусный гепатит, органом-мишенью для которого становится печень. Поэтому окончательно параметры модели подбирались для этого инфекционного заболевания с учётом рекомендаций Вильдановой Ольги Равилевны, старшего лаборанта кафедры микробиологии, вирусологии, иммунологии и клинической лабораторной диагностики ГБОУ ВПО «Южно-Уральский государственный медицинский университет» Минздрава России.
Указанные параметры модели можно рассматривать как иммунологический статус организма. Запишем эти параметры.
Уо — начальная концентрация патогенных антигенов; Го — начальная концентрация антител; С0 — начальная концентрация плазматических клеток; в — коэффициент размножения антигенов; 7 — вероятность нейтрализации антигена антителом; С * — нормальный уровень плазмоклеток; а — вероятность встречи антигена и антитела (эта величина характеризует чувствительность иммунной системы); ^т — обратная величина восстановления органа; ^с — коэффициент времени жизни плаз-моклетки; т — время формирования каскада плазмоклеток (время запаздывания в модели); р — коэффициент производства антител плазмоклеткой; — коэффициент времени распада антитела; п — количество антител, необходимых для распада антигена; М — характеристика здорового органа; М' — характеристика здоровой части органа.
Все коэффициенты модели являются постоянными и положительными.
Приведём теперь несколько теоретических результатов, полученных в работах [1]. Отметим вначале, что система (1) допускает два типа стационарных решений. Одно из них имеет вид
VI = 0, Р = Г*, С1 = С*, шг = 0. (2)
Это стационарное решение описывает состояние здорового организма. Здесь концентрация антигенов и поражённая масса органа равны нулю, а количества плазматических клеток С и антител Г соответствуют иммунологическому статусу здорового человека, т. е. С * и Г *. Показано [1], что достаточным условием асимптотической устойчивости стационарного решения (2) является выполнение неравенства
в < уР*.
Второе стационарное решение имеет вид
V,
Рс(Р/в — 7РС *) в (ар — рс^7) '
г, = в ,
ар/ — п^с72С * 7 (ар — рс^7)
С2
7
^2
т2 =
(3)
Рт
Решение (3) интерпретируется как хроническая форма болезни, если выполняется V, > 0 и С2 > 0.
При компьютерном моделировании [6] особое внимание было уделено случаю, при котором происходит заражение здорового организма малой дозой вирусов В этом случае в начальный момент времени ¿0 отклоняющейся от стационарного решения (2) переменной будет только V(¿), тогда как Г0 = Г*, С0 = С*, т0 = 0. В этой ситуации можно оценить «малость» дозы заражения [7], которая не приводит к потере устойчивости. Такая доза удовлетворяет неравенству
0 < V) < V *
Р/ № * — в)
РТ( '
Величина V* является иммунологическим барьером.
Рис. 1. Результаты моделирования при в < *
На рис. 1 представлены результаты моделирования уравнений системы (1) при в < 7^*. Как показали численные эксперименты, при малых дозах заражения, т. е. дозах заражения, не превышающих иммунологический барьер ("^0 < V*), характер выведения антигенов из организма не зависит ни от дозы заражения, ни от силы
иммунного ответа. Это означает, что такое выведение обеспечивается постоянно присутствующим в организме уровнем антител F *, что соответствует повседневному контакту организма с малыми дозами антигенов. При значительном увеличении дозы заражения по сравнению с величиной иммунологического барьера главную роль в динамике развития заболевания станет играть сила иммунного ответа. Нормальный иммунный ответ способен быстро пресечь развитие инфекции, при слабом иммунном ответе возможна гибель организма. Для приведённого примера иммунологический статус организма задаётся следующим образом: V0 = 0, 005 — начальная концентрация патогенных антигенов; F0 = 120 — начальная концентрация антител; Co = 0 — начальная концентрация плазматических клеток; ß = 15 — коэффициент размножения антигенов; y = 0,4 — вероятность нейтрализации антигена антителом; C* = 10000 — нормальный уровень плазмоклеток; а = 0, 5 — вероятность встречи антигена и антитела (эта величина характеризует чувствительность иммунной системы); ßm = 20 — обратная величина восстановления органа; ßc = 3 — коэффициент времени жизни плазмоклетки; т =14 — время формирования каскада плазмоклеток (время запаздывания в модели); р = 7 — коэффициент производства антител плазмоклеткой; ßf = 5 — коэффициент времени распада антитела; П =10 — количество антител, необходимых для распада антигена; M = 1350 — характеристика здорового органа; M' = 1300 — характеристика здоровой части органа.
0 5000 10000 0 5000 L0000
t t
Рис. 2. Результаты моделирования при ß > y F *
На рис. 2 представлены результаты моделирования уравнений системы (1) при в > yF*. В этом случае протекание заболевания происходит в острой форме, для которой характерно быстрое возрастание количества антигенов в организме до величин, превышающих дозу заражения в несколько раз, а также быстрое выведение антигенов из организма. Такой характер течения болезни обусловлен высоким темпом размножения антигена и сильным, эффективным иммунным ответом.
Для приведённого примера иммунологический статус организма задаётся следующим образом: V0 = 0, 000002 — начальная концентрация патогенных антигенов; F0 = 10 — начальная концентрация антител; C0 = 0 — начальная концентрация плазматических клеток; в = 100 — коэффициент размножения антигенов; Y = 0, 85 — вероятность нейтрализации антигена антителом; C* = 10000 — нормальный уровень плазмоклеток; а = 0, 5 — вероятность встречи антигена и антитела; рт = 20 — обратная величина восстановления органа; рс = 3 — коэффициент времени жизни плазмоклетки; т =14 — время формирования каскада плазмоклеток (время запаздывания в модели); р = 7 — коэффициент производства антител плаз-моклеткой; pf = 5 — коэффициент времени распада антитела; п =10 — количество антител, необходимых для распада антигена; M = 1350 — характеристика здорового органа; M = 1350 — характеристика здоровой части органа.
Список литературы
1. Марчук, Г. И. Математические модели в иммунологии. Вычислительные методы и эксперименты / Г. И. Марчук. — М. : Наука, 1991. — 304 с.
2. Ким, А. В. Математическое и компьютерное моделирование систем с последействием / А. В. Ким, А. Н. Красовский. — Екатеринбург : УГТУ-УПИ, 2010. — 134 с.
3. Марчук, Г. И. Методы вычислительной математики / Г. И. Марчук. — М. : Наука, 1989. — 608 с.
4. Бочаров, Г. А. Численное решение дифференциальных уравнений с запаздывающим аргументом на основе методов Рунге — Кутта — Фельберга / Г. А. Бочаров, А. А. Романюха. — М., 1985. — (Препринт / АН СССР, Отдел вычислит. математики, № 99).
5. Бочаров, Г. А. Численное решение дифференциальных уравнений с запаздывающим аргументом на основе линейных многошаговых методов. Аппроксимация, устойчивость и сходимость / Г. А. Бочаров, А. А. Романюха. — М., 1986. — (Препринт / АН СССР, Отдел вычислит. математики, № 116).
6. Никитина, С. А. Компьютерное моделирование систем с последействием [Электронный ресурс] / С. А. Никитина, А. С. Скорынин. — М. : ФГНУ ИНИПИ РАО. — ОФЭРНиО № 20013 от 24.03.2014.
7. Белых, Л. Н. Качественный анализ простейшей математической модели инфекционного заболевания / Л. Н. Белых, Г. И. Марчук. — Математическое моделирование в иммунологии и медицине. — Новосибирск : Наука, 1982. — С. 5-27.
Поступила в 'редакцию 29.09.2014 После переработки 03.02.2016
Сведения об авторах
Никитина Светлана Анатольевна, кандидат физико-математических наук, доцент кафедры теории управления и оптимизации, Челябинский государственный университет, Челябинск, Россия; e-mail: nikitina@csu.ru.
Скорынин Антон Сергеевич, аспирант математического факультета, Челябинский государственный университет, Челябинск, Россия; e-mail: anton.skorynin@gmail.com.
80
C. A. HuKHTHHa, A. C. CKOPMHHH
Chelyabinsk Physical and Mathematical Journal. 2016. Vol. 1, iss. 1. P. 73-80.
MODELLING OF A DYNAMIC PROCESS WITH AFTER-EFFECT
S. A. Nikitina", A. S. Skoryninb
Chelyabinsk State University, Chelyabinsk, Russia "nikitina@csu.ru, banton.skorynin@gmail.com
Some aspects of the modelling for dynamic systems with delay are considered. Usually it is supposed that the dynamics of a system depends on the current state only. But in the considered systems the behaviour of a process depends on information of its previous states. The arising of the delay in a system essentially changes its dynamics. Special methods are used for modelling of such systems.
Keywords: mathematical modelling, dynamic system with delay, after-effect.
References
1. Marchuk G.I. Matematicheskiye modeli v immunologii. Vychislitel'nye metody i eksperimenty [Mathematical models in immunology. Computational methods and experiments]. Moscow, Nauka Publ., 1991. 304 p. (In Russ.).
2. Kim A.V., Krasovskiy A.N. Matematicheskoye i komp'yuternoye modelirovaniye sistem s posledeystviem [Mathematical and computer modelling of systems with aftereffect]. Ekaterinburg, UGTU-UPI Publ., 2010. 134 p. (In Russ.).
3. Marchuk G.I. Metody vychislitel'noy matematiki [Methods of computational mathematics]. Moscow, Nauka Publ., 1989. 608 p. (In Russ.).
4. Bocharov G.A., Romanyukha A.A. Chislennoye resheniye differentsial'nykh uravneniy s zapazdyvayushchim argumentom na osnove metodov Runge — Kutta — Fel'berga [Computational resolving of differential equations with a retarded argument on the basis of Runge — Kutta — Fehlberg methods]. Moscow, Department of Computational Mathematics of USSR Academy of Sciences Publ., 1985. (In Russ.).
5. Bocharov G.A., Romanyukha A.A. Chislennoye resheniye differentsial'nykh uravneniy s zapazdyvayushchim argumentom na osnove lineynykh mnogoshagovykh metodov. Approksimatsiya, ustoychivost' i skhodimost' [Computational resolving of differential equations with a retarded argument on the basis of linear multi-step methods. Approximation, stability and convergence]. Moscow, Department of Computational Mathematics Academy of Sciences of USSR Publ., 1986. (In Russ.).
6. Nikitina S.A., Skorynin A.S. Elektronnyy resurs «Komp'yuternoye modelirovaniye sistem s posledeystviem» [Electronic resource «Computer modelling of systems with aftereffect»]. Moscow, FGNU INIPI RAO, OFERNiO, no. 20013 from 24.03.2014. (In Russ.).
7. Belykh L.N., Marchuk G.I. Kachestvennyy analiz prosteyshey matematicheskoy modeli infektsionnogo zabolevaniya [Qualitative analysis of a simple mathematical model of infectious disease]. Matematicheskoye modelirovaniye v immunologii i meditsine [Mathematical modelling in immunology and medicine]. Novosibirsk, Nauka Publ., 1982. Pp. 5-27. (In Russ.).
Article received 29.09.2014 Corrections received 03.02.2016
