The research of a mathematical model of immune response reinforcement with application of one-step and multistep methods.pdf Иммунология - это наука об иммунитете живых организмов, а главная задача иммунитета - уничтожение клеток (антигенов), которые генетически отличаются от собственных клеток организма. Вопросами защиты организма от вирусов занимается практическая медицина. Одним из способов борьбы с чужеродными вирусами является стимуляция антителопродуцентами. К настоящему времени в этом направлении накоплен огромный опыт, позволяющий применять в медицине математическое моделирование и численные методы [1-6]. Физико-математическая постановка задачи В настоящей работе исследуется механизм действия САП с применением простейшей математической модели иммунного ответа Г.И. Марчука. Считается, что в организм проникают антигены V(1) (бактерии, вирусы). В процессе выработки антител F(1) участвуют макрофаги M, B-лимфоциты и T-лимфоциты. При попадании в организм антигены поглощаются макрофагами и распознаются T- и B-лимфоцитами. Под действием T-лимфоцитов происходит преобразование B-лимфоцитов в антителопродуценты и формируется популяция антител. Роль антител заключается в связывании антигенов и образовании иммунных комплексов. В данной работе изучаются две возможные гипотезы о механизме действия стимулятора антителопродукции (САП), предложенные в [5]. В рассматриваемом иммунном ответе участвуют антигены V(t), антитела F(t), плазмоклетки C(t) и «молчащие» плазмоклетки B(1). По первой гипотезе, САП вовлекает в антитело-продукцию «молчащие» клетки. Согласно второй гипотезе, после введения САП увеличивается время жизни плазмоклеток. В работе рассматривается математическая модель усиления иммунного ответа [5], которая описывается системой обыкновенных дифференциальных уравнений (ОДУ) с запаздывающим аргументом (задача Коши), в трех вариантах: М0 - контрольная модель (без действия САП), М1 - модель первой гипотезы, М2 - модель второй гипотезы. Уравнения модели М0 с начальными условиями имеют вид dV -= -yFV , dt d- = pC -nrFV F, dC * * * (1) -= aF (t-T )V(t-T )-(c(C-C ), t e(0, T), dt V(0) = V0, С(0) = С*, F(0) = F* = PC-, (f где V(t) - количество неразмножающегося антигена в организме, F(t) - количество антител, C(t) - количество плазмоклеток, производящих антитела. Первое уравнение описывает уменьшение числа антигенов в организме, коэффициент у отражает вероятность нейтрализации антигена антителами. Второе уравнение - производство антител плазмоклетками со скоростью р, уменьшение количества антител после взаимодействия с антигенами с коэффициентом nY и их естественную гибель со скоростью ^f. Третье уравнение характеризует изменение количества плазмоклеток в организме за счет взаимодействия антигенов с антителами с коэффициентом а и старения клеток с константой цс, т является промежутком времени от момента стимуляции лимфоцитов до формирования клона плазматических клеток. Константа C характеризует нормальный уровень плазмоклеток в организме. Обычно начальные условия для уравнений с запаздыванием задаются на интервале [-т ,0]. Однако по биологическому смыслу V(t) = 0 при t < 0, и поэтому начальные условия задаются при t = 0. Предполагается, что все параметры в (1) и V0, C , F неотрицательны. Тогда, согласно [6], решение задачи (1) существует и единственно при всех t > 0. Параметры модели после приведения уравнений к безразмерному виду подбираются из условия устойчивости стационарного решения, о чем речь пойдет ниже. Значение T в безразмерном виде в дальнейшем полагается равным 14. Пусть в соответствии с моделью М1 в организме имеется определенное количество плазмоклеток («молчащие» клетки), которые не могут производить антитела. Для их описания в модель М1 добавляется четвертое уравнение, описывающее изменение количества «молчащих» клеток B(t). Тогда уравнения второй модели с начальными условиями будут иметь вид dV -= -yFV , dt dF — = pC -nyFV -(X F, ^ = ajF (t - t* )V (t -t*) - (c(C - C *) + 5B B, (2) dt dB — = a2F(t -t )V(t -t ) -(BB -SBB, t e(0, T), dt V(0) = V0, C(0) = C*, F(0) = F* =—, B(0) = B0. (f Параметры y, р, nY, ^f, а1, цс в модели М1 взяты такими же, как и в модели М0. В третьем и четвертом уравнениях добавлены члены +5&В и -5в5, описывающие переход «молчащих» клеток в зрелые плазматические клетки. Параметр цв соответствует обратной величине времени жизни B-клеток, а1 Ф а2. В модели М2 до введения САП иммунный ответ развивается в соответствии с моделью М0, поэтому уравнения модели М2 не приводятся. Отличие М2 от М0 состоит в параметре цс - обратной величине среднего времени жизни плазмоклеток, значение которого уменьшается после стимуляции, так как увеличивается время жизни клеток. В моделях М0, Mi, М2 все рассматриваемые величины размерные. На примере модели М0 показывается переход к безразмерному виду. Вводятся обозначения , А * V C F т = uft, Ат = uf т , v = —, c =-, f =-; f f V C F m m m t = —, V = vVm, C = cCm, F = fFm, > m > m > J m > Mf C* Cm = C*, Fm =P—, Vm = —— масштабные множители. Mf ПТ После подстановки (3) в (i) получается система ОДУ dv , Л — = -hf, d т df г Л f =C-f-fv, (4) dc — = И^(т - Ат) f (т - Ат) - h3(c - i), d т v(0) = v0, f(0) = f0, c(0) = C0, C* где h = PY2 , h2 = ap , h3 = — - безразмерные коэффициенты, Дт - безраз-Mf Mf nY Mf мерное время запаздывания, которое в дальнейшем полагается равным шагу сетки. Аналогичным способом приводятся к безразмерному виду уравнения остальных моделей. В случае модели М1 имеем: dv , . d т df г Л "Т = c - f - fv, d т dc — = h2v(т - Ат) f (т - Ат) - h3(c - i) + h4b, (5) d т — = h5v^ - Ат)f (т - Ат) - h6b - h7b, d т v(0) = v„, f (0) = f0, c(0) = ca, b(0) = be, где h2 = ^, h4 = h5 ^, h6 = M h7 Л, ^ = Cm. Mf nY Mf Mf Yn Mf Mf Безразмерный вид модели M2 не приводится, так как она отличается от (4) только параметром h3. (3) Для определения параметров моделей в работе находятся стационарные решения, которые затем исследуются на устойчивость. Стационарные решения выводятся приравниванием правых частей к нулю. Без подробного рассмотрения приводятся стационарные решения изучаемых моделей, полученные в данной работе: vCT = 0, /ст = 1, сст = 1 - для модели М0; vCT = 0, /ст = 1, сст = 1, Ьст = 0 - для модели М1; vCT = 0, /ст = 1, сст = 1 - для модели М2. Исследование устойчивости стационарного решения показывается на примере модели М0. Рассматриваются решения, мало отклоняющиеся от стационарного: 5v = V - Vox , 5/ = / - /ст , 5c = c - Сст . (6) После подстановки (6) в (4) и необходимых преобразований с учетом малости величин второго порядка получается следующая система ОДУ: dr = -К^/ст-/v ст, d5/. = 5c -(1 + VсT)5/ -ист, (7) d т d 5c = ^ 5 / + d т Матрица Якоби системы (7) представима в виде '-К/т -h1vст 0 1 A = /ст -VCT - 1 1 1 h2 /ст h2Vст -h3 J Собственные значения матрицы А на стационарном решении имеют вид Xj = —1; X 2 = —h; Х3 = — h. Все Хг- < 0, i = 1,3, что означает асимптотическую устойчивость по Ляпунову [8] при hi,h3 > 0. Исследование на устойчивость стационарных решений моделей Мь М2 проводится аналогично вышеизложенному. Устойчивость доказана при всех hi > 0, i = 1,7 . При этом показано, что особые точки моделей есть устойчивые узлы. Результаты расчетов Для получения численного решения был проведен вычислительный эксперимент на компьютере. Результаты расчетов оформлены в виде графиков и приводятся по тексту работы. Все рассмотрения проводились для промежутка времени те[0,14] при начальных данных: v0 = 0,04, /0 = 0,065, c0 = 1, hi = 0,5, h2 = 0,45, h3 = 0,2 - для модели М0; v0 = 0,04, /0 = 0,065, c0 = 1, b0 = 0,5, h1 = 0,5, h2 = 0,45, h3 = 0,2, h4 = 0,01, h5 = 0,1, h6 = 0,7, h7 = 0,01 - для модели М1; v0 = 0,04, /0 = 0,065, c0 = 1, h1 = 0,5, h2 = 0,45, h3 = 0,1 - для модели М2. Расчеты проводились с применением одношаговых методов (методы Рунге-Кутты четвертого порядка точности), многошаговых методов (методы Адамса) [7] и чисто неявного А-устойчи-вого метода (метод типа Гира) [7] при n > 40 (расчетный шаг Дт = 14/n). Достоверность программ, по которым велись расчеты, проверялась решением системы ОДУ, имеющей аналитическое решение. Одношаговые методы Рунге-Кутты являются условно устойчивыми, и неустойчивость методов наблюдается при количестве узлов сетки n40 наблюдается сходимость методов. На рис. 1 - 3 представлено поведение антигенов, антител, плазмоклеток и «молчащих» клеток с течением времени в моделях М0, М1, М2 при n = 40. 0 _:_!_ I I I I I 0 2 4 6 8 10 12 т Рис. 1. Поведение v(t), f(т), с(т) с течением времени в модели М^: 1 - v(T); 2 -f(T); 3 - с(т) Рис. 2. Поведение v(t),Дт), с(т), Ь(т) с течением времени в модели Mi: 1 - v(t); 2 - f (т); 3 - с(т); 4 - Ь(т) 0'-;-!- 1 1 1 1 1 0 2 4 6 8 10 12 т Рис. 3. Поведение v(t), f(т), с(т) с течением времени в модели М2: 1 - v(T); 2 - f(T); 3 - с(т) Анализ полученных результатов показывает, что в модели М0 при незначительном изменении количества плазматических клеток количество антител постепенно растет от 0,065 до 1,00509, а количество антигенов убывает (рис. 1). В соответствии с моделью М1 под действием САП «молчащие» клетки созревают в плазматические клетки, вследствие чего количество плазмоклеток увеличивается от 1 до 1,02248 и быстрее растет количество антител от 0,065 до 1,01443 (рис. 2). Согласно второй гипотезе, в модели М2 введение САП приводит к увеличению времени жизни С-клеток. Следовательно, их количество меняется медленнее, чем в модели М0, что способствует соответственно увеличению количества антител (рис. 3). Как следует из графиков, во всех моделях идет вывод антигенов из организма, т.е. выздоровление организма. Однако по сравнению с моделью М0 в моделях М1 и М2 процесс выздоровления организма идет быстрее. Метод Адамса согласно классической теории устойчивости [7] является абсолютно устойчивым методом и не дает ограничений на шаг интегрирования. Процесс вычислений организуется как «предиктор-корректор»: начальные значения предсказываются по методу Рунге-Кутты четвертого порядка точности; затем таблица значений для v(t), f(т), с(т) (модель М0) продолжается по явной и неявной формулам метода Адамса. По явной формуле значения предсказываются, по неявной - уточняются до достижения заданной точности. Использование абсолютно устойчивых методов усложняет процесс вычислений по сравнению с методами Рунге-Кутты, но не требует ограничений на шаг сетки. Численные значения, полученные по методам Адамса и Рунге-Кутты четвертого порядка, мало отличаются между собой и совпадают с точностью до четвертого знака после запятой. Результаты численных расчетов по методу Гира сравнивались с методом Адамса в одноименных точках сетки. Поведение погрешности во всех случаях носит характер «затухающей» волны и максимальное отличие между значениями не превосходит 0,0003. В работе исследовалось действие САП в зависимости от времени его введения. При этом моделировалось введение САП путем изменения параметров h4, h7 в модели М1 и h3 в модели М2 для т = 1, 2, 3, 4, 5. Для удобства дальнейших обсуждений принимается следующий счет безразмерного времени: т = 1 - один день, т = 2 - два дня и т.д. Результаты, представленные на рис. 4 и 5, показывают характер изменения количества плазмоклеток с(т) в зависимости от времени введения САП. Как следует из рисунков, в этом случае модели М1 и М2 предсказывают в рассматриваемом промежутке времени качественно различные динамики течения процесса. В самом деле, в случае модели М1 (рис. 4) в динамике количества зрелых плазмоклеток наблюдается понижение высоты пика в зависимости от времени введения САП: чем позже вводим САП, тем ниже опускается пик. При этом результаты в модели М1 отличаются от результатов, полученных в модели М0 (рис. 6), ростом количества плазмоклеток. В модели М0 максимальное значение с(т) равно 1,020, а в модели М1 оно достигается при введении САП через день (т=1) и равно 1,021. Это пиковое значение появляется на пятый день после иммунизации. Модель М2 дает иную картину динамики С-клеток, и ярко выраженного пика при этом не наблюдается. Однако максимальное значение с(т) равное 1,034 появляется на десятый день и держится два дня, а затем незначительно убывает до 1,033. Результаты качественно похожи между собой и мало отличаются вне зависимости от времени введения фактора (рис. 5). Одной из характеристик роста количества плазмоклеток является величина ЛС, равная разности между количеством плазмоклеток в модели М1 (модели М2) и их количеством в модели М0 (в контроле). Изучалось поведение величины ЛС в моделях М1 и М2 в зависимости от времени введения САП. Изменение данной величины показано на рис. 7, 8. В модели М1 ЛС достигает максимального значения равного 0,0021 на четвертый день при введении САП через день, затем плавно убывает. Аналогичная картина наблюдается и в остальных случаях. Рис. 6. Динамика c(t) в модели Мо Из рис. 7 видно, что более позднее введение фактора приводит к понижению высоты пика величины АС: при введении САП через день она равна 0,0021, через два дня - 0,001 и т.д. Другая картина наблюдается при исследовании величины АС в модели М2 (рис. 8). В данном случае поведение АС не зависит от времени введения фактора и возрастает от 0 до 0,029. Рис. 8. Величина АС для модели М2 В работе исследовалась динамика коэффициента усиления иммунного ответа R(t), равного отношению количества плазмоклеток в соответствии с моделью М1 (моделью М2) к их количеству в модели М0. Характер динамики коэффициента усиления иммунного ответа R(t) с течением времени для модели М1 показан на рис. 9. Из полученных результатов видно, что во всех случаях наблюдается рост значений R(t) от 0 до пикового, затем их плавное убывание. При этом происходит R(t) 1,002 1,0015 1,001 1,0005 1 Рис. 9. Динамика коэффициента усиления иммунного ответа R(t) для модели М1 R(t) 1,025 1,02 1,015 1,01 1,005 1 Рис. 10. Динамика коэффициента усиления иммунного ответа R(t) для модели М2 сдвиг пикового значения вправо в зависимости от времени введения САП и наблюдается его понижение: при введении САП через день максимальное значение равно 1,002 и достигается примерно на четвертый день, при введении САП через два дня оно равно 1,001 и достигается на пятый день и т.д. При введении САП на пятый день значения R(t) приближаются к единице. На рис. 10 изображено поведение коэффициента усиления иммунного ответа R(t) для модели М2. Как видно из рисунка, кривые мало отличаются между собой вне зависимости от времени введения фактора. Значения R(t) меняются от 1 до 1,028, т.е. незначительно отклоняются от единицы. Заключение Таким образом, в работе изучены три математические модели иммунологии с применением численных методов. Из проведенного вычислительного эксперимента на компьютере можно сделать следующие выводы: 1) более раннее введение САП в заболевший организм ускоряет процесс выздоровления; 2) гипотезы, высказанные в работе [5], подтверждены проведенными исследованиями. Из анализа расчетов следует, что введение САП в случае модели М1 дает ярко выраженные результаты на примере поведения АС и R(t). САП в модели М2 почти не зависит от времени ее введения; 3) имеет смысл увеличить интервал наблюдения за характеристиками иммунного ответа, позволяющий увидеть динамику поведения плазмоклеток в случае модели М2.

Скачать электронную версию публикации