Математическая модель влияния иммунотерапии на динамику иммунного ответа mathematical model of immunotherapy effect on the dynamics of immune response Текст научной статьи по специальности «Математика»

У

правление в медицинских и биологических системах

удк 519.622.2

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ВЛИЯНИЯ ИММУНОТЕРАПИИ НА ДИНАМИКУ ИММУННОГО ОТВЕТА

C.B. Русаков, М.В. Чирков

На основе простейшей математической модели инфекционного заболевания построена модель влияния иммунотерапии на динамику иммунного ответа. Сформулирована задача дискретного адаптивного управления, а также предложен алгоритм ее решения. Приведены результаты численного моделирования задачи.

Ключевые слова: математическая модель, адаптивное управление, иммунотерапия.

ВВЕДЕНИЕ

1. ПОСТАНОВКА ЗАДАЧИ

Современный уровень развития иммунологии позволяет рассматривать различные заболевания с единых позиций как процесс взаимодействия иммунной системы с возбудителями болезни. Это дает возможность построения математических моделей абстрактного заболевания, в которых учтены закономерности развития определенного класса болезней. Наиболее существенные характеристики исследуемых процессов отражены в простейшей математической модели инфекционного заболевания, предложенной Г.И. Марчуком [1—3]. Данная модель позволяет изучать влияние внешних воздействий на динамику патологического процесса, обосновать рекомендации по выбору наиболее адекватного лечения. Однако сложность процесса иммунного ответа не позволяет однозначно выбрать критерий управления процессом заболевания непосредственно из содержательных соображений. Поэтому актуальна разработка специализированных методов.

В настоящей работе предпринята попытка развить подход, при котором целью управления может служить обеспечение «идеального» иммунного ответа. Для управления иммунной системой применяется иммунотерапия (введение готовых иммуноглобулинов или донорских антител).

1.1. Простейшая математическая модель инфекционного заболевания

Простейшая математическая модель инфекционного заболевания, предложенная Г.И. Марчуком в 1975 г. [2, 3], описывает фундаментальные механизмы иммунной защиты, сформулированные в клонально-селекционной теории Ф. Берне-та. Модель включает в себя фазовые переменные, которые являются непрерывными функциями:

— концентрацию V(t), част./мл, антигенов (патогенов) в пораженной части органа-мишени;

— концентрацию C(t), клет./мл, плазматических клеток (это популяция носителей и продуцентов антител (иммунокомпетентные клетки);

— концентрацию F(t), част./мл, антител в крови (под антителами понимаются субстраты иммунной системы, нейтрализующие антигены (иммуноглобулины, рецепторы иммунокомпетентных клеток);

— относительную характеристику m(t) пораженного органа.

Модель имеет вид

dV dt

dC dt

= ßV - y FV,

= %(m)aF(t - t) V(t - t) - uXC - C*),

^ = РС - щИГ -МТ = - »тт

м гт

(1)

с начальными условиями при ? е [-т, 0]

Г(?) = г0е(?), С(?) = с*, до = ^ = и*,

»/

(2)

т(?) = 0, где е(?) — функция Хевисайда

е«) = |1 при ' *0 [0 при I< 0.

Биологический смысл параметров модели представлен в табл. 1. Непрерывная невозрастающая неотрицательная функция ^(т) учитывает нарушение нормальной работы иммунной системы вследствие значительного поражения органа. Пусть т* — максимальная доля разрушенных клеток, при которой еще возможна нормальная работа иммунной системы. Тогда функция ^(т) может быть представлена следующим образом:

^(т) =

1, 0 < т < т т - 1

т < т < 1.

(3)

* 1

т - 1

Качественное исследование простейшей модели заболевания позволило выделить четыре основные формы протекания болезни при слабом поражении органа: субклиническую, острую с выздо-

Таблица 1

Параметры простейшей модели инфекционного заболевания

Параметр

Р П

Ъ

а Кт

С

Биологический смысл параметра

Константа скорости размножения антигенов Коэффициент, учитывающий вероятность встречи антигенов с антителами и силу их взаимодействия

Коэффициент стимуляции иммунной системы

Константа скорости естественного старения лимфоцитов

Константа скорости производства антител одной плазмоклеткой

Константа расхода антител на нейтрализацию единицы антигена Константа скорости естественного разрушения антител

Константа скорости разрушения клеток органа-мишени антигеном Константа скорости регенерации органа-мишени

Предсуществующий уровень иммуноком-петентных клеток (плазматических клеток) Время, необходимое для формирования каскада плазмоклеток

ровлением, хроническую и острую с возможным летальным исходом [3].

Остановимся на выборе управления при острой форме инфекционного заболевания.

1.2. Управляемая модель динамики иммунного ответа

Процесс заболевания организма можно рассматривать как управляемый, если понимать под управлением используемые при лечении средства, а также вырабатываемые самим организмом различные гормоны и медиаторы, которые регулируют интенсивности процессов иммунного ответа, а также восстановление пораженных органов и тканей [4].

Одним из эффективных способов лечения острых форм инфекционных заболеваний является иммунотерапия, основанная на введении готовых иммуноглобулинов или донорских антител. В работе [5] предложена модификация простейшей модели инфекционного заболевания с учетом иммунотерапии путем добавления в третье уравнение системы (1) слагаемого, отвечающего за поступление донорских антител.

В предыдущих работах авторов предложен алгоритм дискретного адаптивного управления в простейшей математической модели инфекционного заболевания [6, 7]. В работах [5, 7] собственные и донорские антитела рассматриваются совместно, что, по-видимому, является недостатком.

Будем рассматривать концентрацию собственных и донорских антител отдельно. Для этого введем в модель фазовую переменную, описывающую концентрацию донорских антител в крови К(?), част./мл. Тогда управляемая модель динамики иммунного ответа при инфекционном заболевании будет состоять из пяти обыкновенных дифференциальных уравнений с запаздыванием:

МС

Л

ми ¿к

л

Мт

Л

= рг - уИГ - У1кк,

= ^(т)а(Д? - т) + К(? - т)) Г(? - т)

- »с(С - С*), М = рС - пуИГ - »/Д

= и(/) - ЩУ1*Г - »¿К

= - »тт

(4)

с начальными условиями при ? е [-т, 0]

г(?) = г0е(?), с(?) = с*, до =

к(?) = 0, т(?) = 0

= рС* =

= и *

»/

(5)

а

т

и фазовыми ограничениями

V» 1 0, С(?) 1 0, ДО 1 0, К(?) 1 0, т(?) 1 0.(6)

Функция управления и(?) описывает поступление донорских антител извне и удовлетворяет ограничениям

0 m u(t) m в, 0 m t m t.

(7)

Эти ограничения учитывают физиологически допустимые дозы применения препаратов. Систему уравнений (4) с начальными условиями (5) и ограничениями (6), (7) назовем математической моделью влияния иммунотерапии на динамику иммунного ответа.

1.3. Критерий дискретного адаптивного управления

Сложность процесса иммунного ответа не позволяет однозначно выбрать критерий управления процессом заболевания непосредственно из практических соображений. Трудно надеяться, что когда-нибудь удастся из содержательных соображений обосновать критерий управления, который бы отражал цели управления заболеванием, присущие живому организму [4]. Можно попытаться избежать принципиальных затруднений, связанных с обоснованием критерия управления, если воспользоваться следующим подходом.

В рамках простейшей модели инфекционного заболевания протекание той или иной формы болезни связано с динамикой антигенов V(t). В связи с этим выберем в качестве цели управления обеспечение близости данной характеристики к опорному решению, которое соответствует в некотором смысле «идеальному» иммунному ответу.

Для построения критерия управления рассматриваемым процессом на отрезке [0, T ] зададим равномерную сетку

© = {t; : t. = /At, / = , At = T/N}, (8)

на которой зафиксируем значения концентрации антигенов опорного решения. Для этого рассмотрим задачу минимизации функционала энергетической цены иммунного ответа

E = Lea fF(t - т) V(t - x)dt ^ min , (9)

J a e [ao, ai], т e [0,tq]

где Ь — объем лимфоидной ткани, дренирующей орган-мишень, мл, е — энергетическая цена образования одного лимфоцита, Дж, а0 и т0 — значения соответствующих коэффициентов при естествен-

ном течении заболевания, а а1 — некоторое ограничение сверху на данный коэффициент, так как он не может быть бесконечно большим вследствие биологического смысла задачи. В выражении (9) для каждого а е [а0, ах] и т е [0, т0] значения функций Д? — т) и — т) определяются из решения задачи (1), (2). Пусть минимум функционала Е достигается при некоторых значениях а* е [а0, ах] и т* е [0, т0]. Тогда зафиксируем на сетке (8) множество значений концентрации антигенов, полученное из решения задачи (1), (2) с коэффициентами а* и т*:

V*, / = 1, N.

Будем считать, что условие

V(t.) = V*, / = 1, N

(10)

(11)

соответствует достижению «идеального» иммунного ответа. Таким образом, функция Г(?), описывающая изменение концентрации антигенов, должна проходить через заданный набор точек. В связи с тем, что лечение заболевания может начинаться в разные сроки, а заканчивается при полном выведении антигенов из организма, условие (11) преобразуем следующим образом:

V(tj-) = V;, / е I = {Nx,

N

2>>

1 m n1 m n2 m N.

(12)

Условие (12) будем называть критерием дискретного адаптивного управления процессом иммунного ответа. Таким образом, поставленная задача заключается в построении управления и(?), ? е [0, Т], обеспечивающего выполнение условия (12) при ограничениях (6) и (7).

Будем рассматривать два способа построения управления.

А. Предположим, что донорские антитела вводятся в течение промежутков времени [?,. _ 1, ?,.), / = 1, ..., N с постоянной скоростью, что может соответствовать их поступлению через капельницу. Таким образом, управляющую функцию будем выбирать из множества кусочно-постоянных на отрезке [0, Т] функций

и = {и(?) : и(?) = и,_ 1 е [0, В], ? е [?г _ 1, / = 1, N, и(Т) = им_ 1}.

Б. Дискретное введение донорских антител через интервал времени А?, что может соответствовать инъекциям. В этом случае на каждом отрезке

[t; _ j, t,], i = 1, ..., N, управление входит в начальные условия при t = t. - j следующим образом:

V(t, - !> = к, - 1, c(t, - !> = с, - 1, F(t, - j) = F - 1, K(t, - i) = к - i + u - p m(t, - i) = m - l,

где V - 1, C. - 1, F - 1, K - 1, m, - 1 — значения фазовых переменных на правом конце предыдущего отрезка интегрирования.

Таким образом, в ходе построения дискретного адаптивного управления будем решать последовательность задач Коши на отрезках [t, - 1, t,], i = 1, ..., N, причем значения фазовых переменных в конце предыдущего отрезка служат начальными условиями для следующего.

2. ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧИ АДАПТИВНОГО УПРАВЛЕНИЯ

2.1. Методика расчета

Рассмотрим метод расчета в случае управления А. Будем считать, что при каждом допустимом управлении u = u(t) задача (4), (5) имеет единственное решение

V = V(u(t), t), C = C(u(t), t), F = F(u(t), t), K = K(u(t), t), m = m(u(t), t), u(t) e U,

определенное для всех t e [0, T] и удовлетворяющее условиям (6) и (7).

Сведем решение задачи к последовательному интегрированию уравнений (4) на отрезках [t,. - 1, t,], i = 1, ..., N, на каждом из которых при фиксированном значении u,. - 1 неизвестными функциями являются только фазовые переменные, значения которых на правом конце отрезка служат начальными условиями для следующего промежутка времени.

Для построения управляющей функции u(t) e U применялся следующий алгоритм. На отрезках [t, - 1, t,], i e {1, ..., N}\1, зафиксируем u,. - 1 = 0.

Обозначим ФЦ - 1) = V(u, - 1, t,) — V,*, i e I, где V(u; - 1, t,) — значение функции V(u(t), t) в точке ti при фиксированном значении u, - 1. Множество

значений (10) знается так, что V* < V(0, t,), i e I.

В случае V* ^ V(B, t,) < V(0, t,), i e I, зафиксируем u, - 1 = B. Нахождение управляющей функции на каждом отрезке [t; - 1, tj, i e I, в случае

V(B, t,) < V* < V(0, t,) сводится к решению нелинейного уравнения

Поскольку Ф(0) > 0, а Ф(В) < 0, для решения уравнения (13) применялся метод половинного деления.

Таким образом, управление А можно представить в виде

N -1

и(о = и0 + £ е(? - ^ - 1)(ик - ик - 1). к = 1

В случае управления Б расчет проводится аналогично, только на каждом отрезке вычисляется размер скачка функции К(?), что может описывать инъекции донорских антител через заданные промежутки времени.

2.2. Результаты расчетов — программы лечения

Для численного анализа управления модель заболевания была приведена к безразмерному виду (подобна тому, как это сделано в работе [3]):

^ = — а/ — Л 1 2 9

dS = - т) + - x))v(i - т) - a5(s - 1),

Ф(мг. _ j) = о, о m u - j m в.

(13)

dit = fl4(s - /) - a8/v,

^ = U - a,„kv - a,,k, it 10 11

im

-— = a.v - a-,m, it 6 7

где v = V/ Vm, s = C/C*,f = F/F*, k = K/F*, U = u/F*, Vm — некоторый масштабный множитель для концентрации антигенов, например, биологически допустимая концентрация антигенов в организме (предполагается, что V0 П Vm). В этом случае начальные условия (5) имеют вид

v(t) = Vo9(0, s(t) = 1, f(t) = 1, k(t) = 0, m(t) = 0,

а ограничение (7) задается в виде

о m и(0 m b, о m t m t.

Безразмерные параметры модели описаны в табл. 2. Значения параметров, определяющие острую форму инфекционного заболевания, взяты из работы [3]. В выражении (3) принималось m* = 0,1.

Рассмотрим естественное течение заболевания. Интегрирование задачи (1), (2) проводилось на отрезке времени, длительностью 30 сут. Множество

значений (10) задавалось из решения задачи (1), (2) с параметрами, характеризующими опорное решение, значения которых определяются из выражения (9). В данном случае а3 = 19 000, т = 0, значения остальных параметров не менялись. На рис. 1 представлена зависимость относительной концентрации антигенов от времени. Сплошными линиями показана динамика антигенов при естественном течении заболевания, а штриховыми — в случае опорного решения. Таким образом, управление динамикой иммунной системы направлено на обеспечение иммунного ответа, соответствующего высокой стимуляции иммунной системы при отсутствии запаздывания в реакции на заражение.

Рассмотрим решение задачи дискретного адаптивного управления с шагом А^ = 1 сут. и с шагом А^ = 0,5 сут. в случае описанных способов построения управления. Будем считать, что у{ = у, = п,

= ц-. Зависимость относительной концентрации донорских антител от времени представлена на рис. 2, где по оси ординат отложено отношение концентрации донорских антител к нормальному уровню антител в организме. Сплошными кривыми изображены результаты реализации алгоритма, состоящего в построении дискретного адаптивного управления, в случае, когда сетка (8) задавалась с шагом А/' = 1 сут., а штриховыми — с шагом А/ = 0,5 сут.

В качестве критериев сравнения рассматриваемых программ лечения выберем следующие характеристики:

— степень повреждения органа-мишени, которая определяется как максимальная доля разрушенных антигеном клеток

т

тах т(и(/), /);

111ал I е [ 0, Т]

средняя скорость повреждения организма

Т

I = 1 | аК(и(/), /)Л;

о

Рис. 1. Динамика антигенов

Рис. 2. Концентрация донорских антител: а -

тупление; б — дискретное введение

- непрерывное пос-

Таблица 2

Безразмерные параметры модели

Параметры простейшей модели заболевания Дополнительные параметры управляемой модели

а1 а2 аз а4 а5 аб а7 а8 Т а9 а10 а11 Ь

в * а Ут¥ */С * Н/ Не ^т Нт ПУ^т Уо/Ут У1? * П^т Нк В/¥ *

2 0,8 10 000 0,17 0,5 10 0,12 8 0,5 10-6 0,8 8 0,17 5

Таблица 3

Сравнительные данные программ лечения

Программа лечения mmax I q T

Естественное течение a3 = 104, т = 0,5 2,49 • 10-2 1,03 • 10-3 0 21,59

Опорное решение a3 = 1,9 • 104, т = 0 7,19 • 10-3 3,01- 10-4 0 10,35

Рис. 2, а At = 1 сут 7,27 • 10-3 3,04 • 10-4 3,95 10,45

At = 0,5 сут 6,42 • 10-3 2,72 • 10-4 4,88 9,45

Рис. 2, б At = 1 сут 6,72 • 10-3 2,88 • 10-4 2,60 9,99

At = 0,5 сут 7,17 • 10-3 3,01- 10-4 4,06 10,33

— отношение объема введения донорских антител к нормальному уровню антител в организме, которое в случае управления А определяется как

T

Q = 1 J M(t)dt,

а в случае управления Б — как

N - 1

Q = 1 I

M,,

i = 0

— время полного выздоровления. Будем считать, следуя работе [4], что условие

m(M(t), t) m 5 Vt l T, l 0,

(14)

где 8 > 0 — параметр, выбранный достаточно малым, соответствует практическому выздоровлению с момента Т8, который носит название 8-момента выздоровления.

Количественные характеристики рассмотренных программ лечения приведены в табл. 3, где сначала представлены данные для естественного течения заболевания и опорного решения, а далее — для рассмотренных программ лечения. В условии (14) мы приняли 8 = 0,005. Из сравнения данных характеристик можно сделать вывод о том, что наиболее эффективным способом является дискретное введение донорских антител с интервалом времени, равным одним суткам.

ЗАКЛЮЧЕНИЕ

Результаты расчетов показали, что иммунотерапия позволяет снизить максимальную степень повреждения органа-мишени почти в четыре раза и уменьшить в два раза время, необходимое для полного выздоровления. Рассмотренный способ

управления иммунной системой позволяет достичь иммунного ответа, отвечающего высокой стимуляции иммунной системы при отсутствии запаздывания в реакции на заражение. Предложенный алгоритм также позволяет моделировать различные способы введения препаратов: непрерывное поступление донорских антител и их дискретное введение через определенные промежутки времени.

ЛИТЕРАТУРА

1. Белых Л.Н. Анализ математических моделей в иммунологии / Под ред. Г.И. Марчука. — М.: Наука, 1988. — 192 с.

2. Марчук Г.И. Математические модели в иммунологии. Вычислительные методы и эксперименты. — М.: Наука, 1991. — 304 с.

3. Марчук Г.И. Математические модели в иммунологии. — М.: Наука, 1980. — 264 с.

4. Погожев И.Б. Применение математических моделей заболеваний в клинической практике / Под ред. Г.И. Марчука. — М.: Наука, 1988. — 192 с.

5. Болодурина И.П., Луговскова Ю.П. Оптимальное управление иммунологическими реакциями организма человека // Проблемы управления. — 2009. — № 5. — С. 44—52.

6. Русаков С.В., Чирков М.В. Дискретное адаптивное управление процессом протекания инфекционного заболевания // Вестник Пермского ун-та. Математика. Механика. Информатика. — 2011. — Вып. 1 (5). — С. 84—89.

7. Русаков С.В., Чирков М.В. Дискретное управление в простейшей математической модели инфекционного заболевания // Там же. — 2011. — Вып. 4 (8). — С. 59—63.

Статья представлена к публикации руководителем РРС

В.Ю. Столбовым.

Русаков Сергей Владимирович — д-р физ.-мат. наук,

зав. кафедрой, S (342) 239-65-84, И [email protected],

Чирков Михаил Владимирович — аспирант,

S (342) 242-51-69, И [email protected],

Пермский государственный национальный исследовательский

университет.

0