CC BY
77
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2013 Математика и механика № 3(23)
УДК 519.6
Ш.Х. Султонова, Н.Н. Меркулова
ИССЛЕДОВАНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ УСИЛЕНИЯ ИММУННОГО ОТВЕТА С ПРИМЕНЕНИЕМ ОДНОШАГОВЫХ И МНОГОШАГОВЫХ МЕТОДОВ
Изучаются математические модели усиления иммунного ответа, представляющие собой системы обыкновенных дифференциальных уравнений. Получены стационарные решения в безразмерном виде и показана их устойчивость. Результаты численных расчетов демонстрируют механизм действия стимулятора антителопродукции.
Ключевые слова: математическая модель, стационарные решения, устойчивость, численные методы.
Иммунология - это наука об иммунитете живых организмов, а главная задача иммунитета - уничтожение клеток (антигенов), которые генетически отличаются от собственных клеток организма. Вопросами защиты организма от вирусов занимается практическая медицина. Одним из способов борьбы с чужеродными вирусами является стимуляция антителопродуцентами. К настоящему времени в этом направлении накоплен огромный опыт, позволяющий применять в медицине математическое моделирование и численные методы [1-6].
Физико-математическая постановка задачи
В настоящей работе исследуется механизм действия САП с применением простейшей математической модели иммунного ответа Г.И. Марчука.
Считается, что в организм проникают антигены ¥(() (бактерии, вирусы). В процессе выработки антител _Р(/) участвуют макрофаги М, В-лимфоциты и Т-лимфоциты. При попадании в организм антигены поглощаются макрофагами и распознаются Т- и В-лимфоцитами. Под действием Т-лимфоцитов происходит преобразование В-лимфоцитов в антителопродуценты и формируется популяция антител. Роль антител заключается в связывании антигенов и образовании иммунных комплексов.
В данной работе изучаются две возможные гипотезы о механизме действия стимулятора антителопродукции (САП), предложенные в [5]. В рассматриваемом иммунном ответе участвуют антигены ¥((), антитела _Р(/), плазмоклетки С(0 и «молчащие» плазмоклетки Б(/). По первой гипотезе, САП вовлекает в антитело-продукцию «молчащие» клетки. Согласно второй гипотезе, после введения САП увеличивается время жизни плазмоклеток.
В работе рассматривается математическая модель усиления иммунного ответа [5], которая описывается системой обыкновенных дифференциальных уравнений (ОДУ) с запаздывающим аргументом (задача Коши), в трех вариантах: М0 - контрольная модель (без действия САП), М1 - модель первой гипотезы, М2 - модель второй гипотезы.
Уравнения модели М0 с начальными условиями имеют вид
йУ
----= -уРУ ,
йг
= РС -ПУру р,
йС * * * (1)
— = ар(г -т*)у (г -т*) -цс(С - С ), г е (0, т), йг
у(0) = у0, С(0) = с*, р(0) = р* = рС-,
где У(г) - количество неразмножающегося антигена в организме, р(г) - количество антител, С(г) - количество плазмоклеток, производящих антитела. Первое уравнение описывает уменьшение числа антигенов в организме, коэффициент у отражает вероятность нейтрализации антигена антителами. Второе уравнение - производство антител плазмоклетками со скоростью р, уменьшение количества антител после взаимодействия с антигенами с коэффициентом пт и их естественную гибель со скоростью ц. Третье уравнение характеризует изменение количества плазмоклеток в организме за счет взаимодействия антигенов с антителами с коэффициентом а и старения клеток с константой цс, т является промежутком времени от момента стимуляции лимфоцитов до формирования клона плазматических клеток. Константа С характеризует нормальный уровень плазмоклеток в организме. Обычно начальные условия для уравнений с запаздыванием задаются на интервале [-т ,0]. Однако по биологическому смыслу У(г) = 0 при г < 0, и поэтому начальные условия задаются при г = 0.
Предполагается, что все параметры в (1) и У0, С , р неотрицательны. Тогда, согласно [6], решение задачи (1) существует и единственно при всех г > 0.
Параметры модели после приведения уравнений к безразмерному виду подбираются из условия устойчивости стационарного решения, о чем речь пойдет ниже. Значение Т в безразмерном виде в дальнейшем полагается равным 14.
Пусть в соответствии с моделью М1 в организме имеется определенное количество плазмоклеток («молчащие» клетки), которые не могут производить антитела. Для их описания в модель М1 добавляется четвертое уравнение, описывающее изменение количества «молчащих» клеток В(г). Тогда уравнения второй модели с начальными условиями будут иметь вид
йУ яг/
— = -уРУ ,
йг йр
— = рС -пуру -^ р ,
= ахр(г -т*)У (г -т*) - Цс(С - С*) + 5ВВ, (2)
йг
йВ * *
— = а2р(г -т )У(г -т ) -цВВ -5ВВ, г е (0, Т),
йг
у(0) = у0, С(0) = С*, р(0) = р* =—, В(0) = В0.
Параметры у, р, пу, Ц, а1, цс в модели М1 взяты такими же, как и в модели М0. В третьем и четвертом уравнениях добавлены члены +5в5 и -5в5, описывающие переход «молчащих» клеток в зрелые плазматические клетки. Параметр цв соответствует обратной величине времени жизни В-клеток, а1 Ф а2.
(3)
В модели М2 до введения САП иммунный ответ развивается в соответствии с моделью М0, поэтому уравнения модели М2 не приводятся. Отличие М2 от М0 состоит в параметре цс - обратной величине среднего времени жизни плазмоклеток, значение которого уменьшается после стимуляции, так как увеличивается время жизни клеток.
В моделях М0, Mi, М2 все рассматриваемые величины размерные. На примере модели М0 показывается переход к безразмерному виду.
Вводятся обозначения
, л * V C F
т = -ft, Дт = -f т , v = —, c =------, f =-----;
f f V C F
m m m
t = —, v = vVm, c = cCm, F = fFm,
> ¡fi? m > J m >
-f
C*
Cm = C*, Fm = P—, Vm = —— масштабные множители.
-f ПУ
После подстановки (3) в (i) получается система ОДУ
dv
— = ^^, d т
df г *
Л = c-f-fv, (4)
dc
— = h2v(x - Дт) f (т - Дт) - h3(c - i), d т
v(0) = V0, f(0) = f,, c(0) = c0,
c*
где h = , h2 = ap , h3 = — - безразмерные коэффициенты, Дт - безраз-
-2 Mí пт Mí
мерное время запаздывания, которое в дальнейшем полагается равным шагу сетки.
Аналогичным способом приводятся к безразмерному виду уравнения остальных моделей. В случае модели М1 имеем:
dv , ¿
— = -hlfv, d т
df г *
-f = c - f - fv, d т dc
— = h2 v(т - Дт) f (т - Дт) - h3 (c - i) + h4b, (5) d т
— = h5v^ - Дт)f (т - Дт) - h6b - h7b, d т
v(0) = ve, f (0) = f0, c(0) = ce, b (0) = b0,
где h2 =-aiP-, h4 = ^ h5 =-a2^, h6 = -B, h7 = ^ Bm = Cm.
Mf ПУ -f -f УП -f -f
Безразмерный вид модели M2 не приводится, так как она отличается от (4) только параметром h3.
Для определения параметров моделей в работе находятся стационарные решения, которые затем исследуются на устойчивость.
Стационарные решения выводятся приравниванием правых частей к нулю. Без подробного рассмотрения приводятся стационарные решения изучаемых моделей, полученные в данной работе:
Кт = 0, /ст = 1, сст = 1 - для модели М0;
Кт = 0, /ст = 1, сст = 1, Ьст = 0 - для модели М1;
Кт = 0, /ст = 1, сст = 1 - для модели М2.
Исследование устойчивости стационарного решения показывается на примере модели М0. Рассматриваются решения, мало отклоняющиеся от стационарного:
5у = У - Уст , 5/ = / - /ст , 5с = с - Сст . (6)
После подстановки (6) в (4) и необходимых преобразований с учетом малости ве-
личин второго порядка получается следующая система ОДУ:
‘Г = -/г1^./ст - »15/Vст ,
‘5/. = 5с -(1+ VCT)5/ -ист, (7)
‘ т
‘5с = ЬУст 5 / + Ь5./ст - »,5с.
‘т
Матрица Якоби системы (7) представима в виде
Г-»1./ст -КУст 0 "
4 = -/ст -^т -1 1 .
I »2 /ст »2^ст -»3)
Собственные значения матрицы А на стационарном решении имеют вид
Х1 =-1; X 2 =-^1; Х3 =-»3.
Все Хі < 0, і = 1,3, что означает асимптотическую устойчивость по Ляпунову [8] при Н1,Н3 > 0.
Исследование на устойчивость стационарных решений моделей Мі, М2 проводится аналогично вышеизложенному. Устойчивость доказана при всех » > 0,
і = 1,7 . При этом показано, что особые точки моделей есть устойчивые узлы.
Результаты расчетов
Для получения численного решения был проведен вычислительный эксперимент на компьютере. Результаты расчетов оформлены в виде графиков и приводятся по тексту работы. Все рассмотрения проводились для промежутка времени т е [0,14] при начальных данных: у0 = 0,04, /0 = 0,065, с0 = 1, Н1 = 0,5, Н2 = 0,45, Н3 = 0,2 - для модели М0; у0 = 0,04, /0 = 0,065, с0 = 1, Ь0 = 0,5, И1 = 0,5, Н2 = 0,45, Н3 = 0,2, Н4 = 0,01, к5 = 0,1, Н6 = 0,7, Н7 = 0,01 - для модели Мь у0 = 0,04, /0 = 0,065, с0 = 1, к\ = 0,5, к2 = 0,45, к3 = 0,1 - для модели М2. Расчеты проводились с применением одношаговых методов (методы Рунге-Кутты четвертого порядка точности), многошаговых методов (методы Адамса) [7] и чисто неявного А-устойчи-вого метода (метод типа Гира) [7] при п > 40 (расчетный шаг Дт = 14/п). Досто-
верность программ, по которым велись расчеты, проверялась решением системы ОДУ, имеющей аналитическое решение.
Одношаговые методы Рунге-Кутты являются условно устойчивыми, и неустойчивость методов наблюдается при количестве узлов сетки и<10. Расчеты проводились для различного количества узлов сетки, и было установлено, что при п>40 наблюдается сходимость методов. На рис. 1 - 3 представлено поведение антигенов, антител, плазмоклеток и «молчащих» клеток с течением времени в моделях М0, М1, М2 при п = 40.
Рис. 1. Поведение г(х), f(т), с(х) с течением времени в модели М0 1 - г(т); 2 -Дт); 3 - c(т)
Рис. 2. Поведение v(т),Дт), c(т), Ь(т) с течением времени в модели М1: 1 - v(т); 2 - f (т); 3 - е(т); 4 - Ь(т)
Рис. 3. Поведение v(т),f (т), е(т) с течением времени в модели М2: 1 - v(т); 2 -Дт); 3 - c(т)
Анализ полученных результатов показывает, что в модели М0 при незначительном изменении количества плазматических клеток количество антител постепенно растет от 0,065 до 1,00509, а количество антигенов убывает (рис. 1). В соответствии с моделью М1 под действием САП «молчащие» клетки созревают в плазматические клетки, вследствие чего количество плазмоклеток увеличивается
от 1 до 1,02248 и быстрее растет количество антител от 0,065 до 1,01443 (рис. 2). Согласно второй гипотезе, в модели М2 введение САП приводит к увеличению времени жизни С-клеток. Следовательно, их количество меняется медленнее, чем в модели М0, что способствует соответственно увеличению количества антител (рис. 3). Как следует из графиков, во всех моделях идет вывод антигенов из организма, т.е. выздоровление организма. Однако по сравнению с моделью М0 в моделях М1 и М2 процесс выздоровления организма идет быстрее.
Метод Адамса согласно классической теории устойчивости [7] является абсолютно устойчивым методом и не дает ограничений на шаг интегрирования. Процесс вычислений организуется как «предиктор-корректор»: начальные значения предсказываются по методу Рунге-Кутты четвертого порядка точности; затем таблица значений для v(т), f(т), е(т) (модель М0) продолжается по явной и неявной формулам метода Адамса. По явной формуле значения предсказываются, по неявной - уточняются до достижения заданной точности. Использование абсолютно устойчивых методов усложняет процесс вычислений по сравнению с методами Рунге-Кутты, но не требует ограничений на шаг сетки. Численные значения, полученные по методам Адамса и Рунге-Кутты четвертого порядка, мало отличаются между собой и совпадают с точностью до четвертого знака после запятой. Результаты численных расчетов по методу Гира сравнивались с методом Адамса в одноименных точках сетки. Поведение погрешности во всех случаях носит характер «затухающей» волны и максимальное отличие между значениями не превосходит 0,0003.
В работе исследовалось действие САП в зависимости от времени его введения. При этом моделировалось введение САП путем изменения параметров h4, ^ в модели М! и h3 в модели М2 для т = 1, 2, 3, 4, 5. Для удобства дальнейших обсуждений принимается следующий счет безразмерного времени: т = 1 - один день, т = 2
- два дня и т.д. Результаты, представленные на рис. 4 и 5, показывают характер изменения количества плазмоклеток е(т) в зависимости от времени введения САП. Как следует из рисунков, в этом случае модели М1 и М2 предсказывают в рассматриваемом промежутке времени качественно различные динамики течения процесса.
В самом деле, в случае модели М1 (рис. 4) в динамике количества зрелых плазмоклеток наблюдается понижение высоты пика в зависимости от времени введения САП: чем позже вводим САП, тем ниже опускается пик. При этом результаты в модели М! отличаются от результатов, полученных в модели М0 (рис. 6), ростом количества плазмоклеток. В модели М0 максимальное значение с(т) равно 1,020, а в модели М1 оно достигается при введении САП через день (т=1) и равно 1,021. Это пиковое значение появляется на пятый день после иммунизации. Модель М2 дает иную картину динамики С-клеток, и ярко выраженного пика при этом не наблюдается. Однако максимальное значение е(т) равное 1,034 появляется на десятый день и держится два дня, а затем незначительно убывает до 1,033. Результаты качественно похожи между собой и мало отличаются вне зависимости от времени введения фактора (рис. 5).
Одной из характеристик роста количества плазмоклеток является величина ДС, равная разности между количеством плазмоклеток в модели М1 (модели М2) и их количеством в модели М0 (в контроле). Изучалось поведение величины ЛC в моделях М1 и М2 в зависимости от времени введения САП. Изменение данной величины показано на рис. 7, 8. В модели М1 ЛC достигает максимального значения равного 0,0021 на четвертый день при введении САП через день, затем плавно убывает. Аналогичная картина наблюдается и в остальных случаях.
Рис. 4. Динамика е(т) в зависимости от времени введения САП (модель М^
Ф)
1,03
1,02
1,01
1
0
10
12
Рис. 5. Динамика е(т) в зависимости от времени введения САП (модель М2)
т
Рис. 6. Динамика с(х) в модели М0
Рис. 7. Величина АС для модели М1
Из рис. 7 видно, что более позднее введение фактора приводит к понижению высоты пика величины АС при введении САП через день она равна 0,0021, через два дня - 0,001 и т.д.
Другая картина наблюдается при исследовании величины АC в модели М2 (рис. 8). В данном случае поведение АC не зависит от времени введения фактора и возрастает от 0 до 0,029.
Рис. 8. Величина АС для модели М2
В работе исследовалась динамика коэффициента усиления иммунного ответа Л(х), равного отношению количества плазмоклеток в соответствии с моделью М! (моделью М2) к их количеству в модели М0. Характер динамики коэффициента усиления иммунного ответа R(т) с течением времени для модели М! показан на рис. 9. Из полученных результатов видно, что во всех случаях наблюдается рост значений R(т) от 0 до пикового, затем их плавное убывание. При этом происходит
Рис. 9. Динамика коэффициента усиления иммунного ответа R(т) для модели М!
Рис. 10. Динамика коэффициента усиления иммунного ответа R(т) для модели М2
сдвиг пикового значения вправо в зависимости от времени введения САП и наблюдается его понижение: при введении САП через день максимальное значение равно 1,002 и достигается примерно на четвертый день, при введении САП через два дня оно равно 1,001 и достигается на пятый день и т.д. При введении САП на пятый день значения R(т) приближаются к единице. На рис. 10 изображено поведение коэффициента усиления иммунного ответа R(т) для модели М2. Как видно из рисунка, кривые мало отличаются между собой вне зависимости от времени
введения фактора. Значения R(t) меняются от 1 до 1,028, т.е. незначительно отклоняются от единицы.
Заключение
Таким образом, в работе изучены три математические модели иммунологии с применением численных методов. Из проведенного вычислительного эксперимента на компьютере можно сделать следующие выводы:
1) более раннее введение САП в заболевший организм ускоряет процесс выздоровления;
2) гипотезы, высказанные в работе [5], подтверждены проведенными исследованиями. Из анализа расчетов следует, что введение САП в случае модели М1 дает ярко выраженные результаты на примере поведения АС и R(t). САП в модели М2 почти не зависит от времени ее введения;
3) имеет смысл увеличить интервал наблюдения за характеристиками иммунного ответа, позволяющий увидеть динамику поведения плазмоклеток в случае модели М2.
ЛИТЕРАТУРА
1. Математические модели в иммунологии. Вычислительные методы и эксперименты / ред. Г.И. Марчук. М.: Наука, 1991. 299 с.
2. Математические модели в иммунологии и медицине / ред. Г.И. Марчук. М.: Мир, 1986. 150 с.
3. Марчук Г.И. Математические модели в иммунологии / Г.И. Марчук. М.: Наука, 1985. 240 с.
4. Луговская Ю.П. Математическое моделирование оптимальных процессов лечения инфекционных заболеваний: автореф. дис. ... канд. физ.-мат. наук. Самара, 2009.
5. Математические модели в иммунологии и медицине: сб. статей / пер. с англ. под ред. Г.И. Марчука, Л.Н. Белых. М.: Мир, 1986. 310 с.
6. Белых Л.Н., Марчук Г.И. Качественный анализ простейшей математической модели инфекционного заболевания // Математическое моделирование в иммунологии и медицине. Новосибирскск: Наука, 1982. С. 5-26.
7. Меркулова Н.Н., Михайлов М.Д. Методы приближенных вычислений: учеб. пособие. Томск: ТМЛ-Пресс, 2007. Ч. 2. 288 с.
8. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1968. 720 с.
Статья поступила 20.03.2013 г.
Sultonova Sh.Kh, Merkulova N.N. THE RESEARCH OF A MATHEMATICAL MODEL OF IMMUNE RESPONSE REINFORCEMENT WITH APPLICATION OF ONE-STEP AND MULTISTEP METHODS. Mathematical models of immune response reinforcement as systems of ordinary differential equations are studied. Stationary solutions are obtained in a dimensionless form and their stability is shown. The results of numerical calculations demonstrate the mechanism of the antibody stimulator action.
Keywords: mathematical model, stationary solutions, stability, numerical methods.
SULTONOVA Shahnoza Haidarovna (Tomsk State University)
E-mail: sultonova_sh@sibmail.com
MERKULOVA Nina Nikolaevna (Tomsk State University)
