CC BY
105
УДК 519.622.2
ИДЕНТИФИКАЦИЯ ПАРАМЕТРОВ И УПРАВЛЕНИЕ В МАТЕМАТИЧЕСКИХ МОДЕЛЯХ ИММУННОГО ОТВЕТА
С.В. Русаков, М.В. Чирков
Кафедра прикладной математики и информатики Пермского государственного национального исследовательского университета, Россия, 614990, Пермь, ул. Букирева, 15, e-mail: [email protected]
Аннотация. Цель работы заключается в численном решении задачи дискретного управления иммунным ответом в условиях с неполной информацией. Задача дискретного управления иммунным ответом описана нелинейной системой обыкновенных дифференциальных уравнений с запаздывающим аргументом. Условия с неполной информацией означают, что значения параметров неизвестны, а их оценка корректируется по мере поступления новых экспериментальных данных. Для решения этой задачи предложен алгоритм, позволяющий в рамках математической модели иммунного ответа одновременно строить управление и идентифицировать параметры. Алгоритм строится на основе метода Монте-Карло. Рассматривается управление при острой форме инфекционного заболевания. Идентификация параметров модели проводится по значениям переменных инфекционного процесса. С помощью предложенного алгоритма проведена идентификация параметров базовой математической модели инфекционного заболевания, а также построена программа лечения, основанная на реализации иммунотерапии, которая заключается во введении готовых иммуноглобулинов или донорских антител. Управление выбирается из множества кусочно-постоянных на рассматриваемом отрезке функций. Представленные результаты основаны на имитации экспериментальных значений. В качестве цели управления выбрано обеспечение в некотором смысле «идеального» иммунного ответа, отвечающего высокой стимуляции иммунной системы при отсутствии запаздывания в реакции на заражение. Для этого необходимо динамику антигенов вывести на желаемый режим. Проведено сравнение управления в условиях с неполной информацией с управляющей функцией, построенной с заданным набором значений параметров. Показано, что управление в условиях с неполной информацией лишь незначительно отличается от «идеального» управления. Рассмотренный подход может быть использован для построения различных биомеханических моделей.
Ключевые слова: идентификация параметров, метод Монте-Карло, математическая модель иммунного ответа, дискретное управление, иммунотерапия.
Введение
Математические модели иммунного ответа, широко используемые для исследования динамики иммунной защиты организма при инфекционных заболеваниях [1, 3-7], представляют собой, как правило, нелинейные системы обыкновенных дифференциальных уравнений и содержат большое количество
© Русаков С.В., Чирков М.В., 2014
Русаков Сергей Владимирович, д.ф.-м.н., профессор кафедры прикладной математики и информатики, Пермь
параметров, которые характеризуют иммунный статус организма и свойства антигена. Оценка параметров системы уравнений по клиническим данным позволяет моделировать динамику заболевания у конкретного человека, а также строить прогнозы течения и исхода болезни. Поэтому одна из важнейших задач в этой области заключается в идентификации параметров моделей по данным наблюдений за фазовыми переменными. С помощью оценки параметров можно моделировать управление иммунным ответом и строить программы лечения для конкретного пациента. На практике, как правило, можно измерить значения некоторых фазовых переменных в определенные моменты времени. В связи с этим актуальна разработка эффективных методов идентификации параметров моделей, позволяющих строить управление в процессе идентификации.
Базовая модель инфекционного заболевания
Наиболее общие закономерности функционирования иммунной системы при инфекционных заболеваниях отражены в базовой математической модели инфекционного заболевания, предложенной Марчуком в 1975 г. [5]. Фазовыми переменными модели являются следующие характеристики заболевания: У(() -концентрация антигенов в пораженной части органа-мишени; С(^) - концентрация плазматических клеток; Еф - концентрация антител в крови; ш(() - доля разрушенных антигенами клеток органа-мишени. Базовая модель инфекционного заболевания представляет собой нелинейную систему обыкновенных дифференциальных уравнений с запаздывающим аргументом:
^ = в¥ - уЕУ, сИ
^ = а^(ш) Е ( - т)У ( - т)0(* - т) - ц с (С - С *),
М (1)
СЕ у 4
— = РС - ЩРУ - ц,
сш
~Г = ^ - Ц шш М
с начальными условиями
V(0) = У0, С(0) = Со, Е(0) = Е0, ш(0) = ш0. (2)
Биологический смысл параметров модели представлен ниже. Непрерывная невозрастающая неотрицательная функция ^(ш) учитывает нарушение нормальной работы иммунной системы вследствие значительного поражения органа. Пусть ш* -максимальная доля разрушенных антигенами клеток, при которой еще возможна нормальная работа иммунной системы. Тогда функция ^(ш) может быть представлена следующим образом:
1, 0 < ш < ш*,
4(ш) =
* 1
ш -1
Функция Хевисайда 9(^) определяется по формуле
m -1 . (3)
m < m < 1.
11 при t > 0,
0(t) = 10 t 0 (4)
10 при t < 0.
Параметры базовой модели инфекционного заболевания:
Параметр Биологический смысл
ß Константа скорости размножения антигенов
Y Коэффициент, учитывающий вероятность встречи антигенов с антителами и силу их взаимодействия
а Коэффициент стимуляции иммунной системы
Ц с Константа скорости естественного старения плазматических клеток
Р Константа скорости производства антител одной плазматической клеткой
П Константа расхода антител на нейтрализацию единицы антигена
Ц f Константа скорости естественного разрушения антител
с Константа скорости разрушения клеток органа-мишени антигеном
Ц m Константа скорости регенерации органа-мишени
C * Предсуществующий уровень плазматических клеток
т Время, необходимое для формирования каскада плазматических клеток
Модель (1) описывает общие закономерности, присущие всем инфекционным заболеваниям. Инфекционное заболевание рассматривается как конфликт между патогенным размножающимся возбудителем болезни и иммунной системой организма.
В результате исследования базовой модели инфекционного заболевания выделены качественно отличающиеся друг от друга типы решений, которые интерпретируются как формы протекания болезни (субклиническая, острая, хроническая, летальный исход). Состоянию здорового организма соответствует стационарное решение
V = 0, с1 = С *, Г = = Г*, щ = 0. (5)
^ I
Вид решения однозначно определяется начальными условиями и значениями параметров, которые получили название иммунологического статуса организма. Зная значения параметров модели, можно строить прогнозы течения и исхода болезни для конкретного пациента, а также давать рекомендации по выбору лечения. В связи с этим возникает необходимость в решении задачи идентификации параметров модели по клиническим данным.
Для вычислительных экспериментов модель (1) была приведена к безразмерному виду:
¿V
— = ар -
т
— = а£(ш)/(Х - т)у(Х - т)9(Х - т) - а5($ -1),
аг (6)
— = а4($ - /) -
аХ
¿т
-= а6у - а7 т,
ах 6 7
где V, я, / - относительные концентрации антигенов, плазматических клеток и антител соответственно, V = V / Ут, я = С / С *, / = Е / Е *; Ут - некоторый масштабный множитель для концентрации антигенов, например, биологически допустимая концентрация антигенов в организме.
Моделировалась ситуация заражения здорового организма малой дозой антигенов, поэтому начальные условия задавались в виде
v(0) = V,,, я(0) = 1, /(0) = 1, т(0) = 0. (7)
Для идентификации параметров рассмотрим течение острой формы заболевания, которая характеризуется быстрым ростом концентрации антигенов, выраженным иммунным ответом и резким в силу этого спадом количества возбудителей заболевания до значений, близких к нулю, что понимается как выздоровление.
Значения параметров модели, характеризующие острую форму заболевания, взятые из работы [5], приведены ниже.
Параметр Ü2 a3 a4 a5 a6 a7 a8 T V0 * m
Значение 2 0,8 10 000 0,17 0,5 10 0,12 8 0,5 10-6 0,1
Рис. 1. Динамика переменных инфекционного процесса
f
45 40 35 30 25 20 15 10 5 0
0
10
15
20
25
Рис. 2. Динамика переменных иммунной защиты
30
t, сут
5
На рис. 1 представлена динамика переменных инфекционного процесса при острой форме заболевания. Изменение концентрации антигенов изображено сплошной кривой, а доля разрушенных антигенами клеток органа-мишени -пунктирной. Как видно из рисунка, после полного выведения антигенов из организма требуется большой промежуток времени до полного восстановления органа.
Динамика переменных иммунной защиты при острой форме изображена на рис. 2. Сплошной кривой показано изменение числа плазматических клеток, а пунктирной - концентрация антител. После нейтрализации антигенов еще достаточно продолжительное время наблюдается повышенная концентрация специфичных к данному антигену антител.
Среди приложений математических моделей иммунного ответа важное место занимает задача прогноза течения и исхода болезни для конкретного человека. Применения в клинической практике моделей иммунной защиты организма и подходы к построению управления иммунным ответом представлены в работе [7]. Методы статистического оценивания параметров моделей заболеваний по клиническим данным рассмотрены в работе [4].
Чтобы определить параметры модели при традиционных подходах, необходимо осуществить измерения фазовых переменных в течение всего периода протекания заболевания. Это означает, что полностью оценить параметры можно только к концу заболевания, когда прогноз его течения и исхода теряет свою актуальность. Поэтому необходима разработка таких алгоритмов идентификации параметров, которые бы позволяли осуществлять управление иммунным ответом параллельно с получением экспериментальных значений. Это позволит строить гибкие программы лечения на основе получаемых клинических данных.
Рассмотрим подход, при котором одновременно выполняется управление иммунным ответом и идентификация параметров модели. Будем использовать алгоритм дискретного управления иммунным ответом, предложенный в работах [8, 9]. Идея алгоритма заключается в том, что динамику какой-либо переменной инфекционного процесса необходимо вывести на желаемый режим, соответствующий в некотором смысле «идеальному» иммунному ответу.
Управляемая модель иммунного ответа при инфекционном заболевании может быть представлена в виде
вектор фазовых переменных модели, х(Х) е Яп, п - количество фазовых переменных; управляющая функция и = и(Х), характеризующая скорость введения лекарственных препаратов, выбирается из множества
и = {и (Х): и(Х) = иг-1 е [0, В], Х -1; Х^), ^ = /ДХ, / = 1, Ы, ДХ = , Х0 = 0, и(Т) = иЫ-1}. (9)
Ограничение В > 0 учитывает физиологически допустимые дозы применения препаратов. Параметры модели заданы на множестве
Алгоритм идентификации параметров управляемой модели
(8)
T
A = {а = (а1за2,...,aL): at <ai <а+, i = 1,L}.
(10)
Будем считать, что экспериментальные значения можно получить в узлах сетки
Т
Л = : ^ = iAt, г = г, N М = (11)
где г определяет промежуток времени от момента инфицирования до первого измерения. Соответственно, будем считать и0 = и1 = ... = иг-1 = 0 .
В работе [5] показано, что максимум концентрации антигенов не зависит от дозы заражения, а определяется иммунологическим статусом организма по отношению к данному типу антигена. Поэтому начальную дозу заражения будем считать известной. Также известным будем считать промежуток времени от момента инфицирования до первого измерения.
Алгоритм заключается в следующем. Сначала определяется множество фазовых переменных у(^), по значениям которых будет проводиться идентификация (у (У) е Ят , т < п), и задаются величины допустимого отклонения расчетных значений от экспериментальных данных 8,, , = 1, ..., т. На множестве допустимых значений
параметров А задается сетка дискретизации
а+ - а-
0 = |а : а, =а~ +Д, ] = 0, М,, И =-^т^, г = 1, (12)
Затем случайным образом из узлов сетки (12) задается К наборов значений параметров:
а(к) е 0, к = С~К. (13)
При I еЛ определяются допустимые наборы параметров, т.е. удовлетворяющие следующему критерию идентификации:
у (*,, ^) - УГ ) < 8,, г = г, Ж, к = 1, К, ] = 1, т, (14)
где К1 - количество наборов параметров в момент tг.
Таким образом, после измерения значений фазовых переменных при t = t 1 необходимо найти К{ решений задачи (8) на отрезке [0, ti ], г = г, ..., N и определить допустимые наборы параметров. В качестве оценки параметров при t = t j выбирается среднее значение допустимых наборов:
:Ык) _ _
а;(г) = к=1—, , = 1,ь, г = г, N, (15)
^ г
где 3г - количество допустимых наборов значений параметров в момент времени tj; 3г < Ki, г = г, ..., Ж; К i = 3г-1, г = г + 1, ..., Ж; Кг = К, где К - первоначальное количество наборов параметров; 3 г = К г - , г = г, ..., N где Иг - количество неприемлемых наборов параметров в точке tг.
Обозначим у(а(г), иг, ..., иг, - прогноз значений фазовых переменных на следующий момент времени при данном управлении, где а(г) = {О^-1, а2(г), ..., г = г, ..., N - 1.
В качестве переменной инфекционного процесса, которую необходимо вывести на желаемый режим, выбрана концентрация антигенов, так как с этой характеристикой связано протекание той или иной формы заболевания.
На основе подхода, предложенного в работе [8], определяются соответствующие «идеальному» иммунному ответу значения концентрации антигенов:
X(а(г), 0, ..., 0, I = г, N-1. (16)
Если прогноз уровня концентрации антигенов на следующий момент времени х*(а(|), иг, ..., ы^*, 0, ti+*), I = г, ..., N - 1, не совпадает с «идеальным» значением, то в качестве управления подбирается такая величина и1, I = г, ..., N - 1, которая обеспечивает переход фазовой траектории концентрации антигенов из точки х1эксп (^)
в точку х* (а(|), ti+1), I = г, ..., N - 1.
В качестве решения задачи идентификации параметров берется среднее значение допустимых наборов параметров в конце отрезка интегрирования:
К1 _
а«=, ]=1, (17)
^ N
Таким образом, предложенный алгоритм позволяет построить управление в процессе идентификации параметров модели.
Результаты вычислительных экспериментов
Идентификация параметров модели проводилась по значениям, имитирующим клинико-лабораторные показатели. Для определения значений параметров были выбраны переменные инфекционного процесса: у^) = [у^), т(,)}Т . Множество значений [уэксп ^) = [уэксп ^), тэксп ^)}Т, t еО} задавалось из решения задачи (6), (7) с интервалом времени At = 1 сут. Интегрирование задачи (6), (7) проводилось на отрезке времени Т = 14 сут.
v
t, сут
Рис. 3. Фазовые траектории концентрации антигенов ISSN 1812-5123. Российский журнал биомеханики. 2014. Т. 18, № 2: 259-269 265
Различным значениям параметров модели соответствуют различные фазовые траектории. Идентификация параметров модели проводилась по значениям переменных инфекционного процесса, поэтому фазовые траектории выбранных характеристик должны лежать в некоторой окрестности экспериментальных значений. Таким образом, получаем область, в которой должны находиться фазовые траектории компонентов инфекционного процесса. На рис. 3 для функции у(г) границы этой области изображены пунктирными кривыми. Если при каком-либо наборе параметров траектория выходит за границы этой области в некоторой точке отрезка интегрирования, то дальнейшие вычисления с этим набором параметров не проводятся. На рис. 3 изображены возможные варианты выхода фазовых траекторий за границы допустимой области.
Рассмотрим управление иммунным ответом на примере иммунотерапии, которая заключается во введении готовых иммуноглобулинов или донорских антител. В рамках базовой математической модели инфекционного заболевания реализация иммунотерапии может быть отражена следующим образом [2]:
— = ау - a2Jv,
Ш
— = а3£,(т) / (г - т)у(г - т)9(г - т) - а5 (ь -1), Ш # (18) — = а4(э - /) -а8> + и,
Ш
Шт
-= а6у - а7 т.
Шг 6 7
На рис. 4 изображена реализация предложенного алгоритма. Считалось, что гг = 3 сут. После получения очередного экспериментального значения определяется прогноз уровня концентрации антигенов на следующий момент времени .г1(а(г), иг, ..., и;-1, 0, гг+1), / = г, ..., N - 1. С помощью управления необходимо перевести фазовую траекторию в «идеальную» точку. Как видно из рисунка, фазовая траектория, построенная на основе имитации экспериментальных значений, проходит вблизи «идеальных» точек.
v
г, сут
Рис. 4. Реализация алгоритма: ♦ - «идеальные» значения;--реальная траектория;
▲ - прогноз на следующий момент без управления
Результаты идентификации параметров при К = 104 представлены в табл. 1. Средняя погрешность оценки параметров составляет 8,06%.
Вид управляющей функции представлен на рис. 5. Сплошными линиями показаны результаты построения управления при параллельной идентификации параметров, а штриховыми - управляющая функция, построенная с заданным набором значений параметров. Из рисунка видно, что управление в условиях неполной информации лишь незначительно (до 10%) отличается от «идеального» управления.
В табл. 2 показаны средние значения допустимых наборов параметров в узлах сетки (11). Таким образом, можно сделать вывод о том, что предложенный алгоритм, учитывающий среднее значение допустимых наборов параметров в данный момент времени, позволяет достаточно точно оценить параметры модели и построить управление (программу лечения), близкое к «идеальному».
Таблица 1
Результаты идентификации параметров
Параметр a1 Ü2 a3 a4 a5 a6 a7 a8
Нижняя граница 1,75 0,55 9550 0,145 0,25 7,5 0,095 5,5
Верхняя граница 2,25 1,05 10 550 0,195 0,75 12,5 0,145 10,5
Шаг 0,1 0,1 100 0,01 0,1 1 0,01 1
Оценка 2,150 0,650 10 050 0,175 0,375 10,000 0,112 8,250
Погрешность, % 7,50 18,75 0,50 2,94 25,00 0,00 6,67 3,13
2,0
1,5
1,0
0,5
-Г"
-
1 2 3 4 5 6 7
8 9 t, сут
u
0
Рис. 5. Вид управляющей функции
Таблица 2
Последовательность оценки параметров
ti Hг a1 a2 a3 a4 a5 a6 a7 a8
3 10 000 297 2,000 0,802 10047 0,170 0,503 10,010 0,120 8,006
4 9703 2716 2,000 0,804 10047 0,170 0,502 10,018 0,120 8,008
5 6987 5575 1,993 0,805 10053 0,170 0,505 9,962 0,121 8,008
6 1412 1223 1,983 0,796 10066 0,170 0,493 9,897 0,121 8,093
7 189 105 2,007 0,790 10088 0,170 0,493 9,488 0,121 8,286
8 84 41 2,017 0,799 10115 0,172 0,487 9,616 0,118 8,151
9 43 16 2,061 0,794 10085 0,168 0,472 9,833 0,114 8,093
10 27 8 2,087 0,808 10050 0,167 0,497 10,132 0,112 8,237
11 19 5 2,114 0,800 10014 0,168 0,479 10,143 0,114 8,000
12 14 4 2,120 0,770 9920 0,165 0,480 10,400 0,113 8,100
13 10 3 2,093 0,707 9979 0,172 0,421 10,643 0,109 8,214
14 7 3 2,150 0,650 10050 0,175 0,375 10,000 0,112 8,250
Заключение
Из анализа представленных результатов можно сделать вывод о том, что предложенный алгоритм позволяет достаточно точно оценить параметры модели при относительно небольших вычислительных затратах, а также построить программу лечения в процессе идентификации параметров. Таким образом, данную методику можно применять для идентификации конкретных моделей иммунного ответа при инфекционных заболеваниях по имеющимся клиническим данным. Рассмотренный алгоритм также может быть использован для идентификации других моделей, описанных динамическими системами, при этом необходимо, чтобы были информативны фазовые переменные, по которым проводится идентификация.
Список литературы
1. Белых Л.Н. Анализ математических моделей в иммунологии. - М.: Наука, 1988. - 192 с.
2. Болодурина И.П., Луговскова Ю.П. Оптимальное управление иммунологическими реакциями организма человека // Проблемы управления. - 2009. - № 5. - С. 44-52.
3. Дасгупт Д. Искусственные иммунные системы и их применение. - М.: Физматлит, 2006. - 343 с.
4. Зуев С.М. Статистическое оценивание параметров математических моделей заболеваний. - М.: Наука, 1988. - 172 с.
5. Марчук Г.И. Математические модели в иммунологии. - М.: Наука, 1980. - 264 с.
6. Марчук Г.И. Математические модели в иммунологии. Вычислительные методы и эксперименты. -М.: Наука, 1991. - 304 с.
7. Погожев И.Б. Применение математических моделей заболеваний в клинической практике. - М.: Наука, 1988. - 190 с.
8. Русаков С.В., Чирков М.В. Дискретное управление в простейшей математической модели инфекционного заболевания // Вестник Перм. ун-та. Математика. Механика. Информатика. - 2011. -Вып. 4 (8). - С. 59-63.
9. Русаков С.В., Чирков М.В. Математическая модель влияния иммунотерапии на динамику иммунного ответа // Проблемы управления. - 2012. - № 6. - С. 45-50.
PARAMETER IDENTIFICATION AND CONTROL IN MATHEMATICAL MODELS OF THE IMMUNE RESPONSE
S.V. Rusakov, M.V. Chirkov (Perm, Russia)
The aim of this research is to solve numerically a problem of discrete control of the immune response in conditions with incomplete information. The problem of discrete control of the immune response is described by a nonlinear system of ordinary differential equations with a retarded argument. The conditions with incomplete information imply that the parameter values are unknown and parameter estimation is adjusted as new experimental values become available. In order to solve this problem, we propose the algorithm allowing us to construct a control function and at the same time to identify parameters within mathematical model of the immune response. We have developed the algorithm on the basis of Monte-Carlo method. We deal with the control of immune response in the acute form of infectious disease. By means of the proposed algorithm, we have identified parameters of the basic mathematical model of an infectious disease and constructed the treatment program based on an immunotherapy which consists in introduction of prepared immunoglobulins or donor antibodies. The control is chosen from a set of piece constant functions on considered segment. The presented results are based on the imitation of experimental values. The goal of control is to reach the "ideal" immune response which corresponds to high stimulation of the immune system and the lack of delay in the reaction to infection. For this, it is necessary the antigen dynamics to put into desired level. We have compared the control in conditions of incomplete information with the control function constructed in terms of the given set of parameter values. It is shown that the control in conditions with incomplete information differs only slightly from "ideal" control.
Key words: parameter identification, Monte-Carlo method, mathematical model of immune response, discrete control, immunotherapy.
Получено 8 октября 2013
