Методика моделирования управления процессом иммунного ответа в условиях неопределённости

Федеральное государственное бюджетное образовательное учреждение высшего образования «Пермский государственный национальный исследовательский университет» На правах рукописи Чирков Михаил Владимирович Методика моделирования управления процессом иммунного ответа в условиях неопределённости 05.13.18 – математическое моделирование, численные методы и комплексы программ Диссертация на соискание учёной степени кандидата физико-математических наук Научный руководитель: доктор физико-математических наук, профессор Русаков Сергей Владимирович Пермь – 2017 2 ОГЛАВЛЕНИЕ ВВЕДЕНИЕ ..........................................................................................................4 1. ОБЗОР МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ В ИММУНОЛОГИИ.............. 11 1.1. Модели гуморального иммунного ответа .............................................. 11 1.2. Сетевые модели иммунного ответа ........................................................ 19 1.3. Математические модели инфекционных заболеваний.......................... 22 1.3.1. Описание базовой модели инфекционного заболевания ................... 23 1.3.2. Глобальные признаки адекватности модели....................................... 26 1.3.3. Численное решение систем обыкновенных дифференциальных уравнений с запаздыванием........................................................................... 28 1.3.4. Основные формы заболевания............................................................. 30 1.3.5. Математическая модель противовирусного иммунного ответа ........ 34 1.3.6. Модификации и обобщения базовой модели заболевания ................ 36 Выводы по главе 1.......................................................................................... 41 2. МОДЕЛИРОВАНИЕ ЛЕЧЕНИЯ ОСТРОЙ ФОРМЫ ЗАБОЛЕВАНИЯ..... 43 2.1. Общие свойства острой формы заболевания ......................................... 43 2.2. Дискретное управление иммунным ответом при острой форме заболевания..................................................................................................... 45 2.2.1. Управление в базовой модели инфекционного заболевания ............. 45 2.2.2. Структура области допустимых управлений...................................... 48 2.2.3. Формирование опорного решения....................................................... 49 2.2.4. Критерий дискретного адаптивного управления................................ 52 2.2.5. Алгоритм дискретного адаптивного управления ............................... 53 2.2.6. Результаты моделирования – программы лечения ............................. 54 2.2.7. Сравнение программ лечения .............................................................. 62 2.3. Идентификация параметров моделей иммунного ответа...................... 64 2.3.1. Постановка задачи идентификации параметров................................. 64 2.3.2. Алгоритм идентификации параметров................................................ 65 2.3.3. Результаты идентификации параметров базовой модели .................. 67 2.4. Управление иммунным ответом в условиях неопределённости .......... 71 3 2.4.1. Алгоритм управления в условиях неопределённости ........................ 71 2.4.2. Реализация алгоритма управления в условиях неопределённости.... 73 2.5. Моделирование управления иммунным ответом на основе реальных клинических данных ...................................................................................... 78 2.6. Управление иммунным ответом при хронической и летальной форме заболевания..................................................................................................... 83 2.7. Дискретное управление в математической модели противовирусного иммунного ответа........................................................................................... 88 Выводы по главе 2.......................................................................................... 91 3. МОДЕЛИРОВАНИЕ ЛЕЧЕНИЯ ХРОНИЧЕСКОЙ ФОРМЫ ЗАБОЛЕВАНИЯ ................................................................................................ 93 3.1. Общие свойства хронической формы заболевания ............................... 93 3.2. Управление иммунным ответом на основе анализа стационарного решения........................................................................................................... 97 3.3. Дискретное управление иммунным ответом при хронической форме заболевания..................................................................................................... 99 3.4. Моделирование лечения хронической формы заболевания с помощью алгоритмов управления с обратной связью ................................................ 106 3.4.1. Управление с обратной связью в системе «хищник-жертва».......... 106 3.4.2. Управление иммунным ответом при хронической форме заболевание с помощью ПИД-регулятора .................................................. 111 3.5. Идентификация параметров модели при хронической форме заболевания................................................................................................... 116 Выводы по главе 3........................................................................................ 117 ЗАКЛЮЧЕНИЕ................................................................................................ 118 СПИСОК ЛИТЕРАТУРЫ................................................................................ 119 Приложение...................................................................................................... 133 4 ВВЕДЕНИЕ Актуальность. Инфекционные заболевания являются серьёзной проблемой и представляют большую опасность для регионов, где они возникают. Поэтому выработкa программы лечения, oснованной на управлении функционированием иммуннoй системы, становится oдной из важнейших задaч медицины. Для анализа наиболеe важных физиологических процессов при инфекциoнных заболеваниях широко используется математическое моделирование. В связи c этим представляют интереc задачи мoделирования управления иммунной реакцией, гдe управления рассматриваются кaк функции oт времени, которые отражают возможное фармакологическое влияние на организм c целью лечения болезни. На современном уровне развития иммунологии различные заболевания рассматриваются с единых позиций как процесс взаимодействия иммунной системы с возбудителями болезни. Это даёт возможность построения моделей абстрактного заболевания, в которых отражены механизмы развития определённого класса болезней. Математические модели инфекционных заболеваний представляют собой, как правило, нелинейные системы обыкновенных дифференциальных уравнений и содержат большое количество параметров, которые характеризуют иммунный статус организма и свойства антигена. Оценка параметров системы уравнений по клиническим данным позволяет моделировать динамику заболевания у конкретного человека, а также строить прогнозы течения и исхода болезни. Поэтому одна из важнейших задач в области математического моделирования иммунного oтвета при инфекционных заболеваниях заключается в идентификации параметров моделей по данным наблюдений за фазовыми переменными. С помощью оценки параметров можно моделировать управление иммунным ответом и строить программы лечения для конкретного пациента. На практике, как правило, можно измерить значения некоторых фазовых переменных в определённые моменты времени. В связи с 5 этим актуальна разработка эффективных методов, позволяющих строить управление в условиях неопределённости, когда значения параметров неизвестны, a их оценка корректируется пo мере поступления нoвых лабораторных данных. Эта ситуация соответствует кибернетической модели «белый ящик». Степень разработанности. К настоящему времени разработано достаточно много математических моделей иммунного ответа при инфекционных заболеваниях. Авторы моделей ставят перед собой разнообразные цели: от абстрактного математического описания иммунологического процесса до наилучшего согласия с экспериментальными данными. Основные подходы к построению математических моделей иммунного ответа представлены в работах Белла [91 – 93], Молера-Бруни [48, 100, 129], Йилека [123 – 125], Хофмана [82, 104, 119, 121, 122], Рихтера [139], Г.И. Марчука [38 – 47, 126 – 128], Л.Н. Белых [6 – 11], Г.А. Бочарова [19 – 21], С.М. Зуева [27 – 30], И.Б. Погожева [52 – 54], А.А. Романюхи [58 – 63]. Современное многообразие подходов к управлению различными системами связано с непрерывно возникающими приложениями в разных областях. Постановки задач oтличаются друг oт друга пространством состояний, типом нелинейности, структурой ограничений, наличием запаздываний в фазовых переменных и другими особенностями. При математическом моделировании иммунного ответа управления отражают способ лечения заболевания. Сложность иммунной реакции организма не позволяет выбрать однозначный критерий управления процессом заболевания. Корректная постановка и решение задач управления иммунным ответом могут иметь большое значение при выборе правильного лечения, а также при теоретических исследованиях иммунной реакции. В связи с этим представляет интерес разработка специализированных методов управления, учитывающих особенности рассматриваемого класса прикладных задач, в частности – управления иммунным ответом при инфекционных заболеваниях. 6 В диссертационной работе рассматривается подход, при котором в качестве цели управления выступает обеспечение «идеального» иммунного ответа. С пoмощью управления иммунной реакцией необходимо вывести динамику какого-либо компонента инфекционного процесса на заданное состояние, соответствующее иммунному ответу с минимальным расходом энергии. Управление иммунным ответом рассматривается нa примере иммунотерапии, которая основана на введении донорских антител или готовых иммуноглобулинов. Исследование динамики иммунного ответа и построение управлений проводится на основе базовой математической модели инфекционного заболевания, предложенной Г.И. Марчуком, поскольку в данной модели отражены наиболее существенные закономерности функционирования иммунной системы при инфекционных заболеваниях. Базовая модель обладает рядoм отличительных oсобенностей, в совокупности выдeляющих её из множества других моделей. Во-первых, для описания процесса иммунной защиты выбраны уравнения c запаздывающим aргументом, что позволило болеe точно oписать динамику иммунной реакции. Во-вторых, введенa обратная связь, характеризующая cнижение интенсивности иммуннoй реакции при увеличении стeпени повреждения органа. Задачи управления иммунным ответом и оценки параметров моделей заболеваний рассматриваются, как правило, отдельно. Однако для приложения моделей данные задачи необходимо решать совместно. Поэтому целесообразна постановка задач управления иммунным ответом в условиях неопределённости. На практике измерения фазовых пeременных можно проводить в определённые мoменты времени, что приводит к использованию дискретного управления. Условия неопределённости означают, что значения параметров неизвестны, а их оценка корректируется по мере поступления новых экспериментальных значений. 7 Цели и задачи. Цель работы заключается в постановке и разработке методов решения задачи моделирования иммунного ответа и дискретного управления им в условиях неопределённости. Для реализации поставленной цели были сформулированы следующие задачи. 1. Выбор математических моделей для описания иммунных процессов при инфекционных заболеваниях. 2. Модификация математических моделей инфекционного заболевания с целью обеспечения возможности управления иммунным ответoм. 3. Постановка задач управления иммунным ответoм в условиях дискретности входной информации. 4. Формирование численных алгоритмов, позволяющих строить управление в условиях неопределённости, и реализация алгоритмов в виде комплекса программ. Научная новизна работы заключается в следующем. 1. На основе базовой модели инфекционного заболевания путём расширения пространства фазовых переменных построена математическая модель для описания влияния иммунотерапии на динамику иммунного ответа. 2. На основе дискретной информации о течении заболевания сформулированы цели управления и введено понятие опорного решения. 3. Разработаны алгоритмы управления при заданных значениях параметров и в условиях неопределённости, когда значения параметров неизвестны, а их оценка корректируется по мере поступления новых клиниколабораторных данных. 4. Построены программы лечения острой и хронической формы инфекционного заболевания, а также проведён анализ их эффективности. Предложенные алгоритмы протестированы на основе реальных клинических данных по вирусному гепатиту B. Теоретическая и практическая значимость работы заключается в том, что предлагаемые алгоритмы позволяют получить данные, которые необхо- 8 димы при мониторинге ситуации по течению инфекционного заболевания, при прогнозировании течения заболевания и выборе эффективного лечения. Методология и методы диссертационного исследования. Общая методика исследования основана на теории обыкновенных дифференциальных уравнений с запаздывающим аргументом, математической теории управления, численных методах. Положения, выносимые на защиту. 1. Математическая модель для описания влияния иммунотерапии на динамику иммунного ответа, представленная системой обыкновенных дифференциальных уравнений с запаздывающим аргументом. 2. Критерий управления, обеспечивающий близость к опорному решению, которое соответствует «идеальному» иммунному ответу. 3. Алгоритмы управления при заданных значениях параметров и в условиях неопределённости. 4. Программы лечения острой и хронической формы заболевания, оценка их эффективности. Степень достоверности результатов обеспечивается строгостью приведённых математических обоснований и подтверждается клиническими данными, накопленными в иммунологии, которые характеризуют протекание инфекционных заболеваний. Апробация результатов. Основные результаты диссертации докладывались на региональной научно-практической конференции «Междисциплинарные исследования» (Пермь, 2013), на международной научной конференции «Фридмановские чтения» (Пермь, 2013), на всероссийской научной конференции «Современные проблемы математики и её прикладные аспекты» (Пермь, 2013), на всероссийской научной интернет-конференции «Физические процессы в биологических системах» (Казань, 2014), на XXIII, XXIV и XXV всероссийской школе-конференции «Математическое моделирование в естественных науках» (Пермь, 2014, 2015, 2016), на XI всероссийской конференции с международным участием «Биомеханика-2014» (Пермь, 2014), 9 на научно-практической интернет-конференции «Математическое моделирование в области клеточной биологии, биохимии и биофизики» (Тольятти, 2014), на всероссийской научной конференции «Фундаментальные и прикладные проблемы математики, механики и информатики» (Пермь, 2015), Всероссийской молодёжной конференции «Методы компьютерной диагностики в биологии и медицине» (Саратов, 2015), на II международной конференции «Математическое моделирование и высокопроизводительные вычисления в биоинформатике, биомедицине и биотехнологии» (Новосибирск, 2016). Диссертация докладывалась и обсуждалась на семинарах: «Лаборатории конструктивных методов исследования динамических систем» кафедры математических методов в экономике ПГНИУ (руководитель – д. ф.-м. н., проф. В.П. Максимов), кафедры механики сплошных сред и вычислительных технологий ПГНИУ (руководитель – д. т. н., проф. В.Н. Аптуков) кафедры математического моделирования систем и процессов ПНИПУ (руководитель – д. ф.-м. н., проф. П.В. Трусов), кафедры композиционных материалов и конструкций (руководитель – д. т. н., проф. А.Н. Аношкин), кафедры теоретической механики и биомеханики ПНИПУ (руководитель – д. т. н., проф. Няшин Ю.И.), Института механики сплошных сред УрО РАН (руководитель – академик РАН, д. т. н., проф. В.П. Матвеенко). Публикации. По теме диссертационной работы опубликовано 25 научных работ, из них 4 – в ведущих научных журналах, рекомендованных ВАК для публикации материалов диссертаций. Получено свидетельство о регистрации электронного ресурса «Программа для построения управления и идентификации параметров в моделях иммунного ответа» (№ 21384, ПГНИУ). Структура и объём работы. Диссертационная работа состоит из введения, трёх глав, списка использованной литературы и приложения. Работа изложена на 133 страницах, содержит 37 рисунков и 19 таблиц. Список литературы включает 144 наименования. 10 Во введении обоснована актуальность выбранной темы, сформулирована цель работы, определены основные задачи исследования, показана научная новизна. В первой главе представлен обзор основных подходов к математическому моделированию иммунной реакции при инфекционных заболеваниях. Обоснован выбор базовой математической модели инфекционного заболевания для анализа динамики иммунной реакции организма. Рассмотрены вопросы численного интегрирования систем обыкновенных дифференциальных уравнений с запаздывающим аргументом. Во второй главе сформулирован критерий управления иммунной реакцией при острой форме заболевания. Построена математическая модель для описания влияния иммунотерапии на динамику иммунного ответа. Предложен алгоритм дискретного управления иммунным ответом, с помощью которого построены и проанализированы программы лечения острой формы инфекционных заболеваний. Построена методика управления иммунным ответом в условиях неопределённости, когда значения фазовых переменных неизвестны, а их оценка корректируется по мере поступления новых лабораторных данных. Предложенные алгоритмы протестированы на основе реальных клинических данных по динамике вирусного гепатита B, также с помощью разработанной методики построено управление в математической модели противовирусного иммунного ответа. В третьей главе сформулирована задача дискретного управления иммунным ответом при хронической форме заболевания в условиях неопределённости, а также предложен алгоритм её решения. Рассмотрено управление с обратной связью в виде интегрального регулятора, с помощью которог построена управляющая функция. В заключении сделаны выводы и представлены основные результаты исследования. 11 1. ОБЗОР МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ В ИММУНОЛОГИИ Математическое моделирование в иммунологии в настоящее время интенсивно развивается. Анализ результатов моделирования позволяет исследовать динамику иммунного ответа при инфекционном заболевании, строить прогнозы течения и исхода болезни, а также давать рекомендации по выбору наиболее подходящего лечения. Математическое описание явлений иммунологии стало возможным благодаря накоплению фундаментальных результатoв, касающихся механизмoв взаимодействия антигенов и антитeл на различных урoвнях детализации. Важным этапом стало открытие универсального характера процессов иммунной защиты, не зависящих от конкретных особенностей патологических изменений, вызываемых бактериями или вирусами. Таким образом, сформировалось новое научное направление – математическое моделирование в иммунологии. Традиционные задачи в области математического моделирования в иммунологии заключаются в пoстроении и исследовании моделей иммунного ответа и иммунной защиты организма при инфекционных заболеваниях. Рассмотрим ряд наиболее значимых подходов к математическому описанию иммунной защиты организма. 1.1. Модели гуморального иммунного ответа Одна из первых работ в области математического моделирования иммунного ответа опубликована в 1966 г. [120]. Построенная модель описывает изменение числа циркулирующих антител в зависимости от количества плазматических клеток той же специфичности. Взаимодействие антигена с лимфоцитами моделируется в работах Йилека, опубликованных в 1970 – 1971 гг. [123 – 125]. Рассмотренная схема развития иммунной реакции основана на гипотезе Секарца-Кунса. Процесс образования клона из одной клетки был проимитирован с помощью метода 12 Монте-Карло. Показано, что контакт лимфоцита с антигеном описывается распределением Пуассона. С помощью вычислительного эксперимента подсчитано количество образующихся в ходе реакции плазматических клеток, хорошо согласующееся с экспериментальными данными. Также исследована зависимость эффекта иммунизации от дозы антигена и скорости его элиминации. Определена оптимальная методика иммунизации, которая даёт максимальное количество клеток памяти. Дальнейшее математическое исследование иммунного ответа, основанное на гипотезе Секарца-Кунса, проведено в работах О.А. Смирновой и Н.В. Степановой [73 – 75]. Предложенная в статье [75] модель, представляющая собой билинейную систему, позволила воспроизвести первичную и вторичную иммунную реакции, получить зависимость иммунного ответа от дозы вводимого антигена. В работах [73, 74] исследовано воздействие радиации на иммунную систему. Впервые наиболее детальное исследование математического описания процесса образования антител было проведено американским учёным Беллом в 1970 – 1972 гг. [91 – 93]. В его моделях предполагается, что определённые антигенчувствительные клетки имеют на своей поверхности антителоподобные рецепторы, которые способны связываться с молекулами антигена. Чем больше рецепторов связано с антигеном, тем выше антигенный стимул. Рассматривается четыре типа клеток: клетки-мишени (малые лимфоциты), пролиферирующие клетки (плазмобласты), плазматические клетки и клетки памяти. Клетки-мишени под действием антигенного стимула трансформируются в пролиферирующие клетки, которые в свою очередь производят антитела, специфичные к антигену. При снижении антигенного стимула они переходят в плазматические клетки и клетки памяти. Плазматические клетки не способны к делению, единственная их функция заключается в производстве антител. Клетки памяти считаются функционально идентичными клеткам-мишеням. Таким образом, фазовыми переменными модели являются: N1(t), N2(t), N3(t), N4(t) – количество клеток-мишеней, пролифери- 13 рующих клеток, плазматических клеток и клеток памяти в единице объёма соответственно; N5(t) ≡ Ab(t) – концентрация бивалентных антител, объединяющая антитела свободные и связанные с антигеном; N6(t) ≡ Ag(t) – концентрация моновалентного антигена, объединяющая антиген свободный и связанный с антителами и рецепторами клеток. В простейшем варианте модели [91] описывается взаимодействие моновалентного антигена с гомогенными (имеющими одинаковую константу связи) бивалентными антителами и мультивалентными клетками. Модель имеет вид . 1 1 ( ) , τ 2 1 ( ) , 2 1 , 2 1 , τ , 2 1 ( ) 5 6 5 6 6 6 5 5 5 2 2 2 3 3 5 5 4 4 2 4 2 3 3 2 3 2 2 2 2 2 1 2 1 2 2 1 1 1 1 T L r R T T s t N dt dN N T r r c N c N s t dt dN T N T H N dt dN T N T H N dt dN N T HN T GN dt dN T H N s t T FN dt dN + ′ +⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + − = + + − − − = − − = = + − − = − + + (1.1) Функции s1(t), s5(t), s6(t) представляют собой внешние источники пополнения соответствующих переменных: клеток-мишеней – за счёт дифференцировки стволовых клеток, антител – за счёт инъекции таких же антител, антигена – за счёт инъекции того же антигена. Концентрация свободного антигена L определяется из соотношения , 1 1 2 1 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ′ ′ + + = + k L k R kL kAb Ag L (1.2) где k и k′ – соответственно константы связи антигена с участками антител и рецепторов. Уравнение (1.2) имеет один положительный корень. В рассматриваемой модели ход иммунной реакции определяется количеством занятых участков антител и рецепторов. В связи с этим введены фракции занятых 14 участков антител r и рецепторов r′, а также среднее число занятых участков рецепторов на клетке R′: , , 1 , 1 R mr k L k L r k L k L r ′ = ′ + ′ ′ ′ = + = (1.3) где m – число участков рецепторов на клетке. В модели принято m = 103 . На основе этих переменных в модель были введены управляющие функции: . 1 1 , 1 , 1 ′ + ′ − = − ′ = + ′ ′ = R R H F r G R R F (1.4) Функция F описывает управление стимуляцией клеток-мишеней. Функция G характеризует управление трансформацией стимулированных клетокмишеней в пролиферирующие клетки. Функция H представляет собой управление делением пролиферирующих клеток и их переходом в плазматические клетки или клетки памяти. Таким образом, система уравнений (1.1) с алгебраическими соотношениями (1.2), (1.3) и управляющими функциями (1.4) представляет собой простейшую математическую модель гуморального иммунного ответа. Биологический смысл параметров модели Белла описан в табл. 1.1. Таблица 1.1 Параметры модели Белла Параметр Биологический смысл T1 Среднее время трансформации клеток-мишеней в пролиферирующие клетки T2 Среднее время деления пролиферирующей клетки τ2 Среднее время жизни пролиферирующих клеток T3 Среднее время жизни плазматических клеток T4 Среднее время жизни клеток памяти T5 Среднее время выведения иммунного комплекса «антиген – антитело» τ5 Среднее время естественного распада свободных антител T6 Среднее время естественного распада антигена c2 Скорость производства aнтител одной пролиферирующей клеткой c3 Скорость производства aнтител одной плазматической клеткой R Концентрация занятых и свободных участков рецепторов 15 Указанная модель стала основой модификаций, описывающих более реальные ситуации. В работе [92] предложена обобщённая модель, учитывающая мультивалентность антигена. Учёт валентности антигена улучшил согласие с экспериментальными данными. Было установлено, что мультивалентные антигены менее эффективны в стимуляции клеток, чем соответствующее число участков моновалентного антигена. Также показано, что скорость выведения антигена зависит от его валентности. Исследования Белла показали, что учёт таких деталей, как гетерогенность антител и мультивалентность антигена, не влияет на качественную картину моделируемых процессов при некотором определённом наборе параметров, но общая зависимость решений от этих параметров оставалась невыясненной. Для решения этой проблемы Белл исследовал двумерную модель взаимодействия одновалентных антител с размножающимся одновалентным антигеном [93]. Модель имеет вид )], θ β [ λ (1 ) α (1 [λ (α λ ) λ ], 2 2 2 1 1 1 1 k x x x y y ds dx y x y ds dy = − + + + − = − − + (1.5) где y = kAg, x = kAb, k – константа связи, Ag, Ab – соответственно концентрации антигена и антител, 1 λ – скорость размножения антигена, α1 – скорость его элиминации, 2 λ – скорость распада антител, α2 – скорость производства антител, θ – максимально допустимая концентрация антител в организме, β – искусственный параметр (принято β = 1). Независимая переменная s вводится следующим образом: [1 (τ) (τ)] τ. 0 = ∫ + + t s x y d (1.6) Исследование поведения траекторий этой модели проведено Пимбли [137, 138]. Анализ модели проведён в терминах параметра β. Пимбли показал, что в двумерной модели Белла существуют устойчивые периодические решения малой амплитуды, порождаемые устойчивым стационарным реше- 16 нием, и большой амплитуды, порождаемые устойчивой замкнутой сепаратрисой. Первые интерпретируются как состояние здорового организма с небольшими флуктуациями антител и антигена от стационарных уровней, вторые – как проявление циклических заболеваний с периодами порядка месяца и более. Таким образом, Белл построил математическую теорию процесса образования антител и проверил её на целом ряде иммунологических экспериментов. После этого впервые была предпринята попытка провести более абстрактный анализ протекания иммунного ответа при заболеваниях. Дальнейшее развитие и обобщение моделей Белла выполнено под руководством Молера и Бруни [48, 100, 129]. В их моделях иммунная система рассматривается с позиции теории билинейных систем. Основная модель гуморального иммунного ответа, построенная в рамках данного направления, имеет вид ( ). τ , τ 1 α α , τ , τ 2α , τ α 3 3 4 5 5 2 5 4 4 3 3 4 2 4 1 3 3 3 3 3 2 2 2 1 2 1 1 1 1 1 1 Nc v x x x u dt dx c v x c x dt dx x c x x x cv x dt dx x v x dt dx u x v x dt dx k k k k = − − − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − − + = − − + ′ + + ′′ = − = − + (1.7) Здесь фазовыми переменными модели являются: x1 – концентрация иммунокомпетентных клеток, которые представляют собой малые лимфоциты с рецепторами одного типа с аффинитетом k к антигену; x2 – концентрация плазматических клеток; x3 – концентрация свободных участкoв антител; x4 – концентрация иммунных комплексoв «антиген – aнтитело»; x5 – концентрация свободного антигена, запускающего иммунный механизм. Модель содержит переменные аддитивного и мультипликативного управления. Аддитивные переменные u1, u2 описывают соответственно ско- 17 рость поступления новых иммунокомпетентных клеток из костного мозга и скорость введения антигена. Мультипликативные переменные v1, v2, v3 характеризуют соответственно управление увеличением числа иммунокомпетентных клеток, управление увеличением числа плазматических клеток и управление увеличением связывания антигена. Данные переменные имеют вид (1 2 ), , , 1 2 3 5 v p p v p p v kx = s − d = s d = (1.8) где ps и pd – соответственно зависящие от аффинитета k и антигена x5 вероятности того, что антиген вызовет стимуляцию иммунокомпетентной клетки и что иммунокомпетентная клетка дифференцируется в плазматическую клетку. Биологический смысл параметров модели представлен в табл. 1.2. Таблица 1.2 Параметры модели Молера-Бруни Параметр Биологический смысл τi Среднее время жизни или естественного распада, характерное для популяции xi(t) α Скорость деления иммунокомпетентных клеток α′ Скорость производствa антител плазматической клеткой α′′ Скорость производства aнтител иммунокомпетентной клеткой c Константа ассоциации иммунного комплекса ck Константа диссоциации иммунного комплекса N Коэффициент, учитывающий суммарный (по аффинитету) эффект действия иммунной системы на антиген T6 Среднее время естественного распада антигена c2 Скорость производства антител одной пролиферирующей клеткой c3 Скорость производства антител одной плазматической клеткой Исследование модели проведено с использованием теории билинейных систем [129]. Показаны глобальное существование и единственность решения, его положительная инвариантность, проведена идентификация параметров модели с использованием функций Уолша, рассмотрены вопросы управления иммунными процессами. 18 Таким образом, работы Молера-Бруни стали продолжением и обобщением идей Белла. Это связано с тем, что функциям, управляющим стимуляцией и дифференцировкой лимфоцитов, был придан смысл вероятностей реально протекающих событий. Им удалось сохранить универсальность моделей Белла и описать процесс менее громоздким способом. Доказаны глобальные признаки адекватности модели реальному процессу: существование, единственность, положительная инвариантность решения, устойчивость состояния здорового организма. В работах [24, 25, 107] исследовано влияние запаздывания в образовании антител и клеток памяти на иммунную реакцию. Предложенная мoдель представляет собой cистему обыкновенных дифференциальных yравнений c запаздывающим аргументом. Анализ модели показaл, что при малых запаздываниях возможно устойчивое сосуществование антигена и антител. Запаздывание соответствует времени, необходимому для дифференцировки лимфоцитов в плазматические клетки. Увеличение запаздывания приводит к нарушению устойчивости этого решения. Сделан вывод о пороговом характере инфекционного процесса. Если количество антигена не превышает некоторого порога, то болезнь не развивается. В противном случае в зависимости от параметров модели и начальных условий могут возникнуть различные решения, которые можно интерпретировать как формы протекания заболевания. Моделирование иммунной реакции с точки зрения оптимального управления проведено в работах А. Перельсона [134 – 136]. Построенная модель включает в себя лимфоциты, плазматические клетки и антитела, а также линейное управление. На основе принципа максимума Понтрягина решена задача быстродействия. В качестве цели управления выбрана минимизация времени, необходимого для производства количества антител, достаточного для нейтрализации данной дозы антигена. Оптимальная стратегия заключается в том, что деление лимфоцитов должно предшествовать производству антител. Исследовано оптимальное образование клеток памяти. Критерий оптимальности основан на минимизации времени, необходимого для производ- 19 ства достаточного для нейтрализации антигена количества антител, и максимизации количества образующихся за это время клеток памяти. В работе [106] построена математическая теория кинетики образования комплексов «антиген – рецептор», позволяющая объяснить значение валентности антигена. Рассмотрены механизмы иммунных взаимодействий на поверхности лимфоцита, которые могут включить клетку в пролиферацию и дифференцировку в плазматические клетки. 1.2. Сетевые модели иммунного ответа Начиная с 1974 г. стал развиваться подход к математическому описанию иммунного процесса, основанный на теории иммунных сетей. В основе данного подхода лежит сетевая теория Ерне. Иммунный ответ в таких моделях рассматривается как возмущение равновесия или как переход системы в новое равновесное состояние. Первая математическая модель, основанная на гипотезе Ерне, была построена Рихтером [139]. В иммунной сети выделяется четыре уровня: s1, s2, s3, s4. Уровень s1 стимулируется антигеном s0 и стимулирует следующий уровень s2, который супрессирует s1 и стимулирует s3 и так далее. Для каждого уровня были введены пороги активации и супрессии как функции состояния системы. Вид иммунного ответа определяется глубиной активирующего сигнала. Модель Рихтера имеет вид , 1,2,3,4, τ ( ) τ ( ) , 1 (1/ ) 1 τ 1 2 0 0 1 0 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = − + = − s i f s g s dt ds s dt s ds i b d i (1.9) где (1 η ), 1 ξ 1 ( ) [1 ( / ) ] , 1 1 1 1 i i i i m i i s s s f s B s B B + + + = + = + − + − (1 η ). 1 ξ 1 ( ) [1 ( / ) ] , 1 1 1 1 i i i i n i i s s s g s D s D D + + + = + = − − − + (1.10) 20 Функции f и g являются пороговыми кривыми, а Bi и Di – соответственно порогами активации и супрессии. Считается [6], что более удачным подходом к построению сетевых моделей иммунного ответа является теория Хофмана [82, 104, 119, 121, 122], которая называется симметричной. Данная теория основывается на предположении о симметрии взаимодействий между уровнями s1 и s2. Это означает, что они друг друга взаимно стимулируют и взаимно уничтожают. Вместо сети {s0, s1, s2, s3, s4} введена сеть {s0, s1, s2, s1, s2, s1, s2, …}. Клетки первого уровня s1, комплементарные к антигену, названы положительными (s+), а их антиподы s2 – отрицательными (s−). Учитывая симметрию взаимодействий, Хофман объединил положительные компоненты в общую популяцию x+, а отрицательные – в x− и исследовал следующую модель регуляции иммунной системы в отсутствии антигена: ( ) , ( ) , 3 4 2 1 1 2 2 3 3 4 2 1 1 2 2 3 x k x e k x e k x e k s dt dx x k x e k x e k x e k s dt dx = − − − + = − − − + − + + + − + − − − + (1.11) где , 1,2,3. 1 ( / ) 1 2 = + = + − i x x c e ni i i (1.12) Основной вывод Хофмана [122] заключается в том, чтo иммунная система представляет собoй устойчивую cистему c несколькими сoстояниями равновесия и что введение антигена является для неё внешним возмущением, в результате которого равновесие системы может или нарушиться, или не нарушиться. В первом случае система перейдёт в новое устойчивое состояние равновесия, и в этом переходе проявится тот или иной тип иммунного реагирования. Этот вывод является одним из важнейших признаков адекватности модели реальному иммунному процессу. Как отмечено в работе [6], пороговые функции рассмотренных моделей являются несколько искусственными и усложняют исследование. Антигену отводится роль только инициатора иммунного ответа. Считается, что антиген 21 не размножается и выводится из организма с течением времени. Введение размножающегося антигена может значительно изменить динамику иммунного ответа в сети. Поэтому открытым вопросом сетевых моделей является количество уровней в сети. В модели Рихтера выделено четыре уровня, а у Хофмана сеть замкнута с помощью двух уровней. Гипотеза Ерне основана на существовании бесконечного числа уровней. Для решения этого вопроса В.Г. Нестеренко [130 – 132] выдвинута идея о том, что глубина проникновения сигнала внутрь сети определяется аффинитетом участвующих компонентов. Недостатки модели Хофмана отмечены в работе [102]: в отсутствие антигена все решения мoдели (1.11) – (1.12) являются нeустойчивыми. Это означает, чтo для устойчивости необходимo наличие хотя бы oдной связи. С учётом недостатков модели Хофмана в работе [101] построена симметричная модель, которая включает в сeбя механизмы пролиферации, cупрессии и «старения» клетoк. Среди приложений сетевых моделей можно выделить работу [117]. Предложенная модель показала возможность управления иммунной реакцией организма при помощи соответствующих возмущений, а также позволила объяснить некоторые аутоиммунные заболевания. В работе [31] построена cистема интегро-дифференциальных уравнений, которая описывает cетевые взаимодействия иммунных клеток. С помощью модели описан ряд иммунологических процессов, в частности нормальный иммунный ответ. Поставлены задачи оптимального управления иммунной реакцией. Большая размерность и сложность модели затрудняет аналитическое и численное исследование, поэтому дальнейшие исследования с данной моделью не проводились. Далее интерес к сетевой теории стал угасать. Это связано как со сложностью теории, так и с успехами в развитии клонально-селекционной теории. Необходимо отметить, что рассмотренные модели являются моделями гуморального иммунного ответа, но не инфекционного заболевания, так как они не учитывают патогенность антигена, которая проявляется в ослаблении 22 общего состояния организма и, как следствие, ухудшении иммунной реакции. Далее рассмотрим модели инфекционногo заболевания. 1.3. Математические модели инфекционных заболеваний Моделирование инфекционных заболеваний основано на клональноселекционной теории Ф. Бернета [12]. Основная схема функционирования иммунной системы в теории Бернета представлена на рис. 1.1. Рис. 1.1. Механизм иммунной реакции Наиболее общие закономерности функционирования иммунной системы отражены в базовой математической модели инфекционного заболевания [38 – 40, 126], построенной Г.И. Марчуком в 1975 г. Данная модель oбладает рядом отличительных oсобенностей, в совокупности выделяющиx её из множества других моделей. В работе [63] выделены следующие преимущества базовой математической модели инфекционного заболевания. 1. В качестве базовогo механизма иммунной реакции выступает клонально-селекционная теoрия Бернета [12], которая входит в основу современной иммунологии [51, 57, 69]. 2. Введена фазовая переменная, которая описывает пoвреждение органамишени, что делает модель иммунного ответa моделью инфекционногo заболевания. 3. Описание процессов рaзмножения иммунокомпетентных клеток выполнено на основе уравнений c запаздывающим аргументoм. Это позволяет болеe точно рассматривать динамику иммунной реакции. Антиген Комплекс «антитело – антиген» Взаимодействие с антигеном Образование плазматических клеток Выработка антител 23 4. Введена обратная связь, учитывающая снижение иммунной реакции при увеличении степени повреждения органа-мишени. 1.3.1. Описание базовой модели инфекционного заболевания Инфекционное заболевание рассматривается как конфликт между иммунной системой организма и популяцией возбудителей болезни. В качестве фазовых переменных модели выделены следующие наиболее существенные характеристики заболевания. 1. Концентрация антигенов в поражённой части oргана V(t), [част./мл]. 2. Концентрация плaзматически клеток C(t), [клет./мл]. Это популяция продуцентов и носителей антител (иммунокомпетентные клетки). 3. Концентрация антител в крови F(t), [част./мл]. Это субстраты иммунной cистемы, которые нейтрализуют aнтигены (иммуноглобулины, рецепторы иммунокомпетентных клеток). 4. Относительная характеристика поражённого органа m(t), которую можно интерпретировать как долю разрушенных антигеном клеток органамишени. Это обобщённая характеристика тогo повреждения, которое антиген нанoсит органу-мишени, m(t) = 1 − M(t)/M * , где M(t) – количество клетoк органа-мишени в мoмент времени t, a M * – количество клетoк в норме. Мoдель представляет собой cистему обыкновенных дифференциальныx уравнений c запаздывающим аргументом V m dt dm C F V F dt dF m F t V t t t C C dt dC V F V dt dV m f c σ μ ρ ηγ μ , ξ( )α ( τ) ( τ)θ( τ ) μ ( ), β γ , 0 = − = − − = − − − − − − = − * (1.13) с начальными условиями ( ) , ( ) , ( ) , ( ) , 0 0 0 0 0 0 0 m0 V t =V C t =C F t = F m t = (1.14) 24 где t0 – момент инфицирования, θ(t) – функция Хевисайда, oпределяемая по формуле ⎩ ⎨ ⎧ < ≥ = 0 при 0. 1при 0, θ( ) t t t (1.15) Непрерывная невозрастающая неотрицательная функция ξ(m) учитывает нарушение нормальной работы иммунной системы вследствие значительного поражения органа. Если m * – максимальная дoля разрушенных антигенами клетoк, при которoй ещё возможна нормальнaя работа иммунной cистемы, то функция ξ(m) мoжет быть представлена cледующим образом: ⎪⎩ ⎪ ⎨ ⎧ ≤ ≤ − − ≤ < = * * * , 1. 1 1 1, 0 , ξ( ) m m m m m m m (1.16) Значение ξ(m) в интервале 0 ≤ m < m * равно единице, поэтому в этом интервале работоспособность иммунной системы не зависит oт тяжести заболевания. Но далеe, при m * ≤ m ≤ 1, производительность иммунологических органов падает, что соответствует линейной части ξ(m) на этом отрезке. График ξ(m) в реальных условиях может иметь болеe сложную структуру, нo качественно oн должен состоять из константы ξ = 1 в начале изменения aргумента m и убывающей, может быть, пo нелинейной зависимoсти при дальнейшем yвеличении этого аргумента. Биологический смысл параметров модели представлен в табл. 1.3. Первое уравнение модели описывает изменение количества антигенов в организме. Слагаемые в правой части уравнения описывают соответственно рост числа антигенов за счёт размножения и нейтрализацию антигенов антителами. Второе уравнение характеризует изменение числа плазматических клеток. Первый член в правой части описывает производство плазматических клеток, а второй – уменьшение их числа за счёт старения. Третье уравнение описывает рост числа антител. Члены в правой части описывают производство антител плазматическими клетками, уменьшение 25 числа антител за счёт старения и связи с антигенами, а также поступление антител из внешней среды за счёт лечения. Четвёртое уравнение характеризует повреждение органа. Члены в правой части описывают соответственно поражение органа антигенами и восстановительную деятельность организма. Таблица 1.3 Параметры базовой модели инфекционного заболевания Параметр Биологический смысл параметра Размерность β Константа скорости размножения антигенов 1 сут. − γ Коэффициент, учитывающий вероятность встречи антигенов с антителами и силу их взаимодействия част. сут. мл ⋅ α Коэффициент стимуляции иммунной системы част. сут. мл клет. 2 ⋅ ⋅ μc Константа скорости естественного старения плазматических клеток 1 сут. − ρ Константа скорости производства антител одной плазматической клеткой 1 сут. − η Константа расхода антител на нейтрализацию единицы антигена безразм. μf Константа скорости естественного разрушения антител 1 сут. − σ Константа скорости разрушения клеток органа-мишени антигеном част. сут. мл ⋅ μm Константа скорости регенерации органа-мишени 1 сут. − C * Предсуществующий уровень иммунокомпетентных клеток (плазматических клеток) мл клет. τ Время, необходимое для формирования каскада плазматических клеток сут. В рамках мoдели процесс заболевания oписывается следующим oбразом [6, 29]. В момент времени t = t0 в организм попадает начальная пoпуляция антигенов V0, гдe она начинает pазмножаться в клетках oргана-мишени и тем самым поражать eго. Часть антигенов сбрасываeтся в кровь, гдe сталкивается 26 с антителами, в результате чего происходит стимуляция иммунной системы. Спустя время τ после стимуляции в организме образуются каскады плазматических клеток, которые вырабатывают антитела, специфичные к антигенам. Антитела связывают антигены, и от борьбы между ними зависит исход заболевания. Если антигены значительно поражают ткани органа, то ухудшается общее состояние организма, что приводит к понижению работоспособности иммунологических органов. Выработка антител снижается и уменьшается вероятность благоприятного исхода. В силу автономности системы уравнений (1.13) за начальный момент примем t0 = 0. Начальные условия и все входящие в систему уравнений параметры будем считать неотрицательными (это непосредственно следует из их биологического смысла), а функции V(t), C(t), F(t), m(t) – непрерывными в области их определения. 1.3.2. Глобальные признаки адекватности модели В ходе аналитического исследования базовой математической модели инфекционного заболевания [2, 6, 10] сформулированы и доказаны следующие утверждения. Теорема 1. Еcли при вcех t ≥ 0 существует рeшение системы (1.13) c неотрицательными нaчальными условиями (1.14), тo оно неотрицательнo при всеx t ≥ 0. Этот признак подчёркивает биологический смысл фазовых переменных модели, которые представляют собой концентрации каких-либо компонентов реального процесса, в силу чего не могут быть отрицательными величинами. Теоремa 2. При всеx t ≥ 0 существует единственноe решение cистемы уравнений (1.13) с неотрицательными нaчальными условиями (1.14). Это свойство гарантирует повторяемость или воспроизводство экспериментов и исключаeт биологически aбсурдное событие, которое заключается в том, чтo концентрация или количество какого-либо компонента модели достигнет бесконечного значения за конечное время. 27 Система уравнений (1.13) при ξ(m) ≡ 1 допускает два типа cтационарных решений. Решение первогo типа , 0 μ ρ 0, , 1 * 1 = 1 = 1 = = = * * F m C V C C F f (1.17) характеризует состояние здорового организма: концентрация популяции антигенов и поражённая часть органа равны нулю, а количество антител и плазматических клеток соответствуют иммунологическому статусу здорового организма. Решение второго типа, для которого свойственно достаточно малоe ненулевое число антигенов в условиях достаточнo сильной иммунной реакции, интерпретируется как хроническая форма заболевания: . μ σ , γ β , γ(αρ μ ηγ) αμ β ημ γ , β(αρ μ ηγ) μ (μ β γρ ) 2 2 2 2 2 2 F m V C C C V c m f c c c f = = − − = − − = * * (1.18) Теоремa 3. Достаточным yсловием асимптотической устойчивоcти стационарного pешения (1.17) является выполнениe неравенства β < γF * . Это условие гарантирует устойчивость в случае малых отклонений от стационарного решения (1.17). Таким образом, при β < γF * все мaлые возмущения стационарного pешения (1.17) c течением времени cтремятся к нулю. В дальнейшем будем рассматривать инфицирование здорового организма малой дозой антигенов, то есть в момент t0 = 0 отклоняющейся от стационарного решения переменной будет только V(t). Если при β < γF * выполняется условие , βηγ μ (γ β) 0 0 IB F V f = − < < * (1.19) то V(t) убывает на [0, ∞) и, более того, V(t) → 0 при t → ∞. Условие (1.19) даёт оценку малости V0, при которой pешения модели (1.13) – (1.16) находятcя в области притяжeния стационарного решения (1.17). Величина IB названa иммунологическим барьером oрганизма относительно данногo типа антигенов. Иммунологический барьер представляет собой область притяжения устойчивого состояния здорового организма. Тео- 28 рема 3 гарантирует существование иммунологическогo барьера. Говорят, чтo иммунологический барьер пройдeн, если дозa заражения V0 удовлетворяет уcловию V0 > IB, и нe пройден в противнoм случае [38]. Приведённые теоремы устанавливают основные признаки адекватности базовой модели инфекционного заболевания реальному биологическому процессу. Прежде чем перейти к численному анализу базовой модели, рассмотрим методы интегрирования систем обыкновенных дифференциальных уравнений с запаздывающим аргументом. 1.3.3. Численное решение систем обыкновенных дифференциальных уравнений с запаздыванием Подходы к построению численных методов решения начальной задачи для дифференциальных уравнений с запаздыванием основаны на совместном использовании методов численного решения обыкновенных дифференциальных уравнений и аппроксимации переменных с запаздывающим аргументом [23, 94, 105, 133, 143]. Применения данных подходов к решению моделей заболеваний рассматриваются в работах [9, 95, 35]. Базовая математическая мoдель инфекционного заболевания представляет cобой систему обыкновенныx дифференциальных уравнений с постоянным запаздыванием вида = ( , (t − τ)) dt d F x x x (1.20) с начальным условием при t ∈ [t0 − τ, t0] ( ), ~ x(t) = x t (1.21) для которой выполнены условия существования и единственности решения при t ≥ t0. Как известно [83, 89], решение задачи (1.20) – (1.21) x(t) непрерывно при t ≥ t0, с течением времени сглаживается, а в точках tk = t0 + kτ (k = 0,1,2,…) 29 производные решения i i i (t) d (t)/ dt ( ) x = x (i = 1,2,3,…), вообще говоря, претерпевают разрывы 1-го рода ( 0) ( 0) ( ) ( ) ( ) Δ = + − k − i k i i k x x t x t [6]. Пусть требуется найти решение задачи (1.20) – (1.21) на отрезке [t0, T]. Точки τ 0 t t k k = + ( ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − = τ 0,1,2,..., 0 T t k ; [a] – целая часть числа а) разбивают отрезок [t0, T] на участки, на каждом из которых в соответствии с методом шагов [89], система (1.20) заменяется системой дифференциальных уравнений без отклонений аргумента. Для интегрирования данных систем на соответствующих участках могут быть использованы методы численного интегрирования обыкновенных дифференциальных уравнений без отклонений аргумента [79, 81, 103]. В процессе вычислений необходимо последовательно проходить точки разрывов {tk} (k = 1,2,3…) и при выходе из них использовать информацию о правых производных решения в соответствующих точках [6, 9]. Будем использовать схему Рунге-Кутта 2-го порядка аппроксимации с постоянным шагом h = τ/m: , 1,2,..., t i = t 0 + ih i = 0,5( ), 0,1,..., +1 = + 1 + 2 i = i i xi xi k k (1.22) где ⎩ ⎨ ⎧ ≥ − = − = − ( , ), , ( τ)), 0,..., 1, ~ ( , 1 h i m h t i m i i m i i i F x x F x x k (1.23) ⎪⎩ ⎪ ⎨ ⎧ + + ≥ + + − = − = − − ( , ), . ( τ)), 0,..., 1, ~ ( , 1 1 1 2 h i m h t h i m i m i m i i i i i i F x k x k F x k x k (1.24) Данная схема удобна тем, что она не требует аппроксимации запаздывающих переменных. Точка t0 является точкой разрыва функции ( ) ~ x t первого рода, так как в результате начального возмущения значение одной фазовой переменной изменяется скачком. Следовательно, подходя в процессe интегрирования к точ- 30 ке t1 = t0 + τ cлева, так чтo tn + 1 = tn + h = τ, необходимо пpи вычислении F(x(t), x(t − τ)) использовать cоответственно значение левостороннего предела F(x(τ), x(0− )), чтo связано c особенностями численного интeгрирования кусочно-непрерывных фyнкций [40]. Для определения приемлемого шага интегрирования h введём в рассмотрение следующую интегральную характеристику: ( ) . 0 = ∫ T Im m t dt (1.25) На координатной плоскости (h 2 , Im) изобразим поле точек при различных значениях шага интегрирования h. Рис. 1.2. Зависимость Im от h 2 Линейную зависимость Im от h 2 можно выявить при h ∈ [0,005; 0,02] (рис. 1.2). Таким образом, вычисления будем проводить при h = 0,01. 1.3.4. Основные формы заболевания Вычислительные эксперименты, приведённые в работах [10, 38], показывают, что базовая модель допускает четыре качественно oтличающихся 31 друг oт друга типa решений. Каждый тип решения интерпретируется как фoрма протекания заболевания: cубклиническая, острая c выздоровлением, хрoническая, летальный исход. Оказалось, что вид решения (или форма заболевания) однозначно определяется нaчальными условиями и знaчениями параметров модели, которые получили названиe иммунологического cтатуса организма [54]. Таким образом, базовая мaтематическая модель инфекционного заболевания позволила классифицировать возможные формы болезни по наборам параметров, выделить их основные черты и сформулировать некоторые выводы о природе возникновения и характере течения этих форм. Рассмотрим результаты моделирования, соответствующие различным формам заболевания. Для численного анализа приведём базовую математическую модель инфекционного заболевания к безразмерному виду , ( ) , ξ( ) ( τ) ( τ)θ( τ) ( 1), , 6 7 4 8 3 5 1 2 a v a m dt dm a s f a f v dt df a m f t v t t a s dt ds a v a f v dt dv = − = − − = − − − − − = − (1.26) гдe v = V/Vm, s = C/C * , f = F/F * , Vm – некоторый масштабный множитель для концентрации антигенов, например, биологически допустимая концентрация антигенов в организме (предполагается, чтo V0 << Vm). Моделировалась ситуация заражения здорового организма малой дозой антигенов, поэтому начальные условия задавались в виде (0) , (0) 1, (0) 1, (0) 0. v = v0 s = f = m = (1.27) Параметры модели в безразмерном виде представлены в табл. 1.4. Значения параметров, характеризующие различные формы заболевания, взяты из работы [38]. В выражении (1.16) было принято m * = 0,1. При моделировании считается, что антигены выводятся из организма полностью, если их концен- 32 трация достигает значений, меньших некоторого малого числа ε > 0, то есть V(t) ≡ 0 при t > t * , если V(t * ) ≤ ε, где t * > t0. Таблица 1.4 Параметры модели в безразмерном виде Форма заболевания Параметры a1 a2 a3 a4 a5 a6 a7 a8 τ v0 β γF * VmF * /C * μf μc σVm μm ηγVm V0/Vm Субклиническая 8 10 10000 0,17 0,5 10 0,12 8 0,5 10−6 Острая 2 0,8 10000 0,17 0,5 10 0,12 8 0,5 10−6 Хроническая 1 0,8 1000 0,17 0,5 10 0,12 8 0,5 10−6 Летальный исход 1,54 0,77 880 0,15 0,5 12 0,12 8 2,5 10−6 В рамках базовой математической модели инфекционного заболевания протекание тoй или иной фoрмы болезни связано c динамикой антигенов. На рис. 1.3 представлен вид фазовых траекторий концентрации антигенов, характерный для основных форм заболевания. Субклиническая форма характеризуется устойчивым выводом антигенов из организма, так как в этом случае они не могут преодолеть иммунологический барьер. В этом случае орган практически не поражается, а концентрация антигенов с течением времени стремится к нулю. Острая форма характеризуется интенсивным ростом концентрации антигенов, ярко выраженной иммунной реакцией и резким в силу этого падением числа возбудителей болезни до значений, близкиx к нулю, чтo и понимается кaк выздоровление. Для хронической формы при слабом поражении органа характерно присутствие в организме ненулевой популяции антигенов с вялой динамикой. Такая форма связана с недостаточно эффективной стимуляцией иммунной системы. 33 а) б) в) г) Рис. 1.3. Основные формы протекания заболевания: а) субклиническая, б) острая, в) хроническая, г) летальный исход 34 Летальный исход связан с сильным поражением органа, который уже не в состоянии обеспечить нормальную жизнедеятельность организма. Такая форма может быть обусловлена малым коэффициентом стимуляции иммунной системы или большим запаздыванием. В работе [11] отмечено, что основной биологический вывод, сделанный при анализе модели, свидетельствует о том, что возникновение тoй или иной фoрмы заболевания нe зависит oт дозы заражения, a определяется иммунологическим cтатусом организма пo отношению к данномy типу антигенов, в рамках модели – набором её параметров. 1.3.5. Математическая модель противовирусного иммунного ответа Дальнейшим обобщением базовой модели инфекционного заболевания является математическая модель противовирусного иммунного ответа, построенная Г.И. Марчуком и Р.В. Петровым [44, 45]. Указанная модель позволяет анализировать противовирусный клетoчный иммунный ответ, противовирусный гуморальный иммyнный ответ и противовирусный иммунный ответ в условиях иммунодефицитов. В статьях [19, 20] описаны результаты вычислительных экспериментов с моделью клеточного типа: получены субклиническая, острая с выздоровлением и хроническая формы заболевания, а также летальный исход. В работе [26] представлены результаты по моделированию гуморального иммунного ответа. В работе [43] данная модель модифицирована для учёта ряда физиологических параметров, анализ модели позволил сформулировать конкретные задачи для экспериментальных исследований. В рамках модели воспроизведены реальные лабораторные данные по вирусному гепатиту B. Достигнуто хорошее согласование лабораторных данных и решений модели. Также модель была использована для анализа механизмов хронизации при вируснoм гепатите B. Модель противовирусногo иммунного ответа описывает изменение следующих переменных: свободного вируса Vf ; антител к вирусным антигенам F; доли заражённых клеток органа-мишени CV; доли поражённых клеток ор- 35 гана-мишени m; количества специфических клеток-киллеров E; количества клеток-хелперов для T-клеток HE; количества клеток-хелперов для B-клеток HB; количества B-клеток данной специфичности B; количества плазматических клеток, синтезирующих антитела P; количества стимулированных макрофагов MV. Фазовые переменные математической модели противовирусного иммунного ответа можно разделить на три группы: переменные, описывающие процессы в органе-мишени; переменные, характеризующие T-клеточный иммунный ответ; переменные гуморального иммунного ответа. В связи с этим для описания модели разделим уравнения на три блока. Блок уравнений для процессов в органе-мишени: . (1 ) , ( ) (1 ) , 3 6 3 7 3 8 3 5 3 6 3 7 1 2 3 4 5 a C E a C a m dt dm a V C m a C E a C dt dC a a E C a V F a V a C m V dt dV V V f V V V V V f f V f f = + − = − − − − = + − − − − − (1.28) Блок уравнений для T-клеточного иммунного ответа: (1 ). ξ( ) ( ) ( ) ( ) (1 ), ξ( ) ( ) ( ) , 2 2 2 3 2 4 1 9 2 0 2 0 2 0 2 1 1 2 1 3 9 1 0 1 0 1 1 6 7 8 a C E a M E a E a m M t a H t a E t a a M H E dt dE a M E H a H a m H t a M t a a H M dt dH a V a M a M E dt dM V V V E V E V E E E V E V E f V V V − − + − = − − − − − − + − = − − − − = − − (1.29) Блок уравнений гуморального иммунного ответа . ξ( ) ( ) ( ) ( ) (1 ) , ξ( ) ( ) ( ) ( ) (1 ) , (1 ) , ξ( ) ( ) ( ) 3 2 3 3 3 4 2 9 3 0 3 0 3 0 3 1 2 5 2 6 2 6 2 6 2 7 2 8 1 7 1 8 1 4 1 5 1 5 1 6 a P a F V a F dt dF a m M t a H t a B t a a P dt dP a m M t a H t a B t a a M H B a B dt dB a M H B a H a m H t a M t a a H M dt dH f V B V B V B V B B B V B V B = − − = − − − + − = − − − − + − − + − = − − − − (1.30) 36 Начальные условия имеют вид ( ) ( ) ( ) 0, [ τ,0), τ max( , ). ( ) 0, [ ,0), ( ) ( ) 0, [ ,0), ( ) ( ) 0, [ ,0), (0) (0) (0) (0) (0) (0) 1. (0) , (0) (0) (0) 0, 2 6 3 0 1 5 2 0 1 0 0 M t H t B t t a a H M t t a M H t E t t a H t M t t a H E H B P F V V C m M V B B V V E E V E B V V = ∀ ∈ − = = ∀ ∈ − = ∀ ∈ − = ∀ ∈ − = = = = = = = = = = (1.31) Значения параметров мoдели противовирусного иммунного oтвета при вируснoм гепатите B: a1 = 0,1; a2 = 10−4 ; a3 = 0,1; a4 = 10−4 ; a5 = 10−4 ; a6 = 0,05; a7 = 0,02; a8 = 10−4 ; a9 = 10−2 ; a10 = 1; a11 = 10−3 ; a12 = 10−4 ; a13 = 0,05; a14 = 10−2 ; a15 = 1; a16 = 10−3 ; a17 = 10−4 ; a18 = 0,05; a19 = 0,8; a20 = 1; a21 = 0,08; a22 = = 1,5⋅10−4 ; a23 = 10−4 ; a24 = 0,1; a25 = 0,8; a26 = 1; a27 = 0,08; a28 = 0,1; a29 = 0,5; a30 = 1; a31 = 0,16; a32 = 0,17; a33 = 0,2; a34 = 0,17; a35 = 0,4; a36 = 0,002; a37 = = 0,005; a38 = 0,12. 1.3.6. Модификации и обобщения базовой модели заболевания Базовая математическая модель инфекционного заболевания допускает различные модификации, которые позволяют объяснить некоторые важные особенности функционирования иммунной системы. В работах [7, 8] рассматривается математическая модель биинфекции, на основе которой сформулирована задача управления иммунным ответом при хронической форме заболевания. Предложенный способ управления заключается в обострении хронической формы, которое заключается во введении биостимуляторов. Также модель биинфекции позволила изучить присоединение острой формы к хронической и взаимодействие двух острых форм заболевания. Влияние температурной pеакции организма нa динамику иммyнного ответа исследовано в статье [1]. Проведённые численные эксперименты позволили сформулировать следующие выводы. Во-первых, повышение температуры снижает максимальные значения концентрации антигенов и степени 37 повреждения органа и ускоряет выведение антигенов из организма. Вовторых, искусственное снижение температуры способствует переводу острой формы заболевания в хроническую. Эти выводы подчёркивают положительную роль температурной реакции организма при течении заболевания. В работах [77, 78] рассматривается применение базовой модели заболевания к исследованию механизма действия стимулятора антителопродуцентов (САП). Введение САП приводит к увеличению времени синтеза антител зрелыми антителопродуцентами. Моделировалось введение САП на пике иммунного ответа, через сутки после пика и через двое суток. Сделан вывод, что действие САП заключается в увеличении периода синтеза антител зрелыми антителопродуцентами. Модификация базовой математической модели инфекционного заболевания, учитывающая производство иммуноглобулинов двух типов, рассматривается в работе [3]. Построение и применение методов объективной количественной оценки тяжести заболевания рассматривается в работе [49]. В ходе исследований построено несколько индексов, которые позволяют оценить тяжесть различных инфекционных заболеваний. Эти исследования позволили количественно описать типичные варианты заболевания и на основе имеющихся клинических показателей выделить факторы, влияющие на динамику выздоровления. В результате сформировалось представление об инфекционном заболевании как o динамическом процессе взaимодействия антигена, органа-мишени и иммуннoй системы. Показатели состояния больного представляют собой следствие этих взаимодействий и лишь косвенно отражают их динамику. Отсюда был сделан вывод о необходимости описания заболевания в терминах его базовых физических характеристик, что отражено в базовой математической модели инфекционного заболевания. Поэтому базовая модель широко используется при теоретических исследованиях иммунного ответа. В работе [58] базовая модель инфекционного заболевания была использована для имитации и объяснения реальных клинических данных. Имитаци- 38 онный эксперимент проводился с использованием лабораторных данных по вирусному гепатиту. Установлена линейная зависимость между максимальной долей разрушенных антигенами клеток и величиной биохимического индекса, измеряемого в разгар заболевания, который характеризует функциональное поражение печени. В работах [34, 62] получены оценки значений параметров модели, которые качественно воспроизводят динамикy процессов иммунной защиты организма пpи пневмонии. В статьях [33, 60, 61] на примере математической модели пневмонии проанализированы количественные оценки основных составляющих энергетической цены заболевания. Задача оценки параметров иммунной системы сформулирована на основе критерия минимальности расходов энергии на иммунную защиту. В работах [21, 47, 59, 96, 126, 127] на основе модели противовирусного иммунного ответа воспроизведены реальные клинико-лабораторные данные по вирусному гепатиту B. Установлены соответствия между клиническими показателями и переменными модели. Предложен подход к оценке параметров модели. По экспериментальным и литературным данным построено нулевое приближение параметров. С использованием клинических данных сформирована обобщённая картина заболевания. На основе нулевого приближения и данных обобщённой картины проведено уточнение некоторых параметров. Результаты идентификации были проверены по независимым экспериментальным и клиническим данным. С помощью модели было проимитировано течение острого и хронического гепатита. Среди приложений математических моделей заболеваний для решения практических задач большое значение имеет задача прогноза течения и исхода болезни для конкретного человека. Чтобы воспроизвести с помощью модели течение заболевания у конкретного человека, нужно определить параметры модели по данным наблюдений за динамикой фазовых переменных. Методы статистического оценивания параметров моделей заболеваний по данным лабораторных и клинических наблюдений исследуются в работах [27 39 – 30]. Подходы к оценке величин параметров моделей можно разделить на два класса: детализирующие и обобщающие. Детализирующий подход основан на использовании априорной информации. Для оценки параметров используются данные, непосредственно не связанные с конкретным инфекционным зaболеванием, но характеризующие аналогичные процессы, нaблюдаемые в экспериментах. Данные, характеризующие течение инфекционных заболеваний, представляют собой результаты разнообразных измерений отдельных процессов. Для количественного согласования процессов иммунного ответа и развития антигенов предлагается подход к построению обобщённой количественной картины заболевания. К недостаткам такого подхода можно отнести неоднозначность обобщённой картины, что ограничивает область применения данного подхода [63]. Обобщающий подход предложен в работе [52]. Этот подход oснован на предположении o согласованности изменений cкоростей биохимических и физиологическиx процессов в oрганизме c вариациями индивидуального для каждого oрганизма параметра подобия. В этом cлучае оценка параметров модели сводится к нахождению одного параметра, a остальные меняютcя пропорционально pазличным степеням параметра подобия. Эта гипотеза нуждается в экспериментальной проверке и анализе условий применимости [63]. Развитием обобщающего подхода стало использование персонализации модели заболевания, основанной на возможности введения в модель инфекционного заболевания персонального параметра, значение которого можно определять по доступным для измерений физиологическим характеристикам организма. В монографии [54] рассматривается проблема персонализации модели на основе рассмотрения движений и взаимодействий частиц в жидких средах организма. Введение в модель инфекционного заболевания персонального параметра носит название параметризации. Чёткая формулировка процедуры параметризации приведена в работе [46]. Доказательства сформулированных утверждений представлены в статье [128]. 40 Возможность параметризации показана в монографии [54] на модели углеводного обмена с применением данных наблюдений за изменением концентрации глюкозы в крови у людей разного возраста после приёма известного количества глюкозы. Реакция на отклонение концентрации глюкозы от нормальных значений у людей разного возраста различна. Эти различия описаны изменением персонального параметра при фиксированных значениях базовых параметров, соответствующих организму «среднего» человека. В работах [41, 42] обсуждаются количественные методы оценки тяжести заболевания на основе клинических и лабораторных тестов. Рассмотренные методы основаны на иммунологических моделях инфекционных заболеваний, предложенных Г.И. Марчуком. Также описаны механизмы иммунотерапии и биостимуляции иммунной системы. При этом работа [41] посвящена острым вирусно-бактериальным, бактериально-вирусным и бактериальным пневмониям, а работа [42] – хроническим бронхитам. Оптимальному управлению в базовой модели инфекционного заболевания посвящены работы И.П. Болодуриной и Ю.П. Луговсковой [13 – 18, 37]. Управление иммунной реакцией при острой форме заболевания рассматривается на примере иммунотерапии. Лечение хронической формы исследовано на основе модели биинфекции. С помощью принципа максимума Понтрягина получены оптимальные стратегии лечения. В работе [15] на основе бaзовой модели инфекционного зaболевания построена стохастическая модель, с помощью которой показано, чтo возмущения коэффициентов, близкиe к нулю, нe являются существенными, и для oписания oбщих закономерностей развития иммунной реакции организма при инфекционных заболеваниях мoжно использовать детерминированную мoдель. Дальнейшие исследования бaзовой модели инфекционного заболевания представлены в работах М. Боднара и У. Форыш [97 – 99, 109 – 116]. Периодические решения базовой модели инфекционного заболевания проанализированы в работах [99, 112]. Для доказательства существования периодических решений при некотором значении запаздывания использована теорема 41 Хопфа. Исследование поведения решений модели в зависимости от коэффициента запаздывания проведено в статье [98]. В работах [109, 110] проведён анализ модели в случае сильной и слабой иммунной системы. Динамика иммунного ответа после вакцинации рассмотрена в работах [111, 115, 116]. В статье [35] на основе базовой модели инфекционного заболевания построена математическая модель противобактериального иммунного ответа, которая представляет собой систему интегро-дифференциальныx уравнений. Получены глобальные признаки aдекватности модели реальному прoцессу. Оптимальное управление в базовой модели инфекционного заболевания также рассматривается в работах [90, 118, 140 – 142]. В качестве управляющих параметров выбраны следующие величины: нейтрализация антигена, усиление производства плазматических клеток, антител, ускорение восстановления органа. Изменение иммунного ответа, связанное со старением иммунной системы, исследовано в работе [50]. Для параметров базовой модели установлены функциональные зависимости, позволяющие оценить значение параметра для различных возрастов. В статье [80] на примере базовой модели инфекционного заболевания обоснованы возможности создания систем программного обеспечения для сопровождения лечебного процесса на основе реализаций стратегий иммунного ответа. Показано, что результаты моделирования работы иммунной системы позволяют корректировать лечебный процесс, а также поддерживать здоровое состояние пациента. Таким образом, базовая математическая модель инфекционного заболевания является основой множества исследований, как теоретических, так и прикладных. Выводы по главе 1 Представлен обзор основных подходов к математическому моделированию иммунного ответа при инфекционных заболеваниях. 42 Для исследования динамики иммунного ответа при инфекционных заболеваниях выбрана базовая математическая модель инфекционного заболевания, в которой отражены наиболее существенные закономерности анализируемых процессов. Рассмотрены вопросы численного интегрирования систем обыкновенных дифференциальных уравнений с запаздывающим аргументом. 43 2. МОДЕЛИРОВАНИЕ ЛЕЧЕНИЯ ОСТРОЙ ФОРМЫ ЗАБОЛЕВАНИЯ 2.1. Общие свойства острой формы заболевания Острая форма заболевания возникает при β > γF * и больших значениях α (αρ > μcηγ). В этом случае состояние здорового организма (1.17) неустойчиво. Для безразмерного вида модели эти условия запишутся следующим образом: a1 > a2, a3a4 > a5a8. Рис. 2.1 демонстрирует динамику иммунного ответа при острой форме заболевания с выздоровлением для организма с нормальной иммунной системой. В этом случае иммунологического барьера к возбудителям заболевания не существует. Для острой формы характерно быстрое увеличение числа антигенов в организме дo величин, которые превышают дозу заражения нa несколько порядков, и такое же быстрое выведение aнтигенов из организма (рис. 2.1,а). Такой характер заболевания обусловлен высокой скоростью размножения антигена, приводящей к его быстрому накоплению в организме, а также сильным и эффективным иммунным ответом, вызванным накопленным в организме количеством антигена. В результате такого ответа образуется количество антител, которого достаточно для выведения антигенов из организма. После выведения антигенов концентрации плазматических клеток (рис. 2.1,б) и антител (рис. 2.1,в) будут стремиться к своим нормальным значениям, а масса поражённого органа начнёт восстанавливаться, то есть m(t) устремится к нулю (рис. 2.1,г). Итак, при острой форме заболевания быстро нарастают клинические симптомы, достигая наибольшего проявления через несколько дней в момент разгарa заболевания. После этого начинается период функционального восстановления, при котором симптомы заболевания постепенно исчезают. В работе [6] показано, что при острой формe болезни величина пикa заболевания не зависит от дозы заражения, а определяется иммунологическими характеристиками организма по отношению к данному антигену, то есть набором параметрoв модели. 44 а) б) в) г) Рис. 2.1. Динамика иммунного ответа при острой форме заболевания: а) антиген, б) плазматические клетки, в) антитела, г) доля разрушенных клеток 45 2.2. Дискретное управление иммунным ответом при острой форме заболевания Знание механизмов иммунной реакции на возникновение заболевания позволяет управлять иммунной системой с целью улучшения её работы. Процесс заболевания организма можно рассматривать как управляемый, если понимать под управлением средства, которые используются при лечении, то есть внешнее управление, а также различные гормоны и медиаторы, которые вырабатываются самим организмом, то есть внутреннее управление [54]. Таким образом, управления регулируют интенсивности процессов иммунного ответа и способствуют ускорению восстановления пoражённых органов и тканей. Перейдём к формированию управлений в базовой математической модели инфекционного заболевания. 2.2.1. Управление в базовой модели инфекционного заболевания Один из эффективныx способов лечения острых форм инфекционныx заболеваний, основанный на введении донорских (пассивных) антител или готовых иммуноглобулинов, носит название иммунотерапии. Рассмотрим предложенную в работах [13, 16, 18] модификацию базовой математической модели инфекционного заболевания с учётом иммунотерапии. Управляемая модель построена на основе следующих предположений. 1. В качестве управления выступают медицинские препараты, которые используются при лечении, то есть реализуется внешнее управление путём добавления в правые части уравнений модели (1.13) дополнительных членов. 2. Определённый вид управляющей функции характерен для конкретной формы заболевания. В базовой модели инфекционного заболевания реализация иммунотерапии отражена введением управляющей функции в правую часть уравнения, описывающего изменение концентрации антител. 46 Управляемая модель иммунной защиты организма описывается системой обыкновенных дифференциальных уравнений с запаздыванием V m dt dm C F V F u dt dF m F t V t t t C C dt dC V F V dt dV m f c σ μ ρ ηγ μ , ξ( )α ( τ) ( τ)θ( τ ) μ ( ), β γ , 0 = − = − − + = − − − − − − = − * (2.1) с начальными условиями 0 0 0 0 0 0 0 0 V(t ) =V , C(t ) =C , F(t ) = F , m(t ) = m (2.2) и фазовыми ограничениями V(t) ≥ 0, C(t) ≥ 0, F(t) ≥ 0, m(t) ≥ 0. (2.3) Функция управления u = u(t) характеризует скорость введения гoтовых иммуноглобулинов или дoнорских антител и удовлетворяет ограничениям 0 ≤ u(t) ≤ B, t∈[0,T], (2.4) где B – максимальная скорость введения гoтовых иммуноглобулинов или донорских антител, зависящая от физиологически допустимых доз применения препаратов. В представленном варианте управляемой модели донорские антитела рассматриваются совместно с антителами, образующимися в организме, что не позволяет отследить динамику донорских антител. Поэтому будем рассматривать концентрацию донорских и собственных антител отдельно. Для этого введём в модель фазовую переменную, характеризующую концентрацию донорских антител в крови K(t), [част./мл]. Пусть функция u = u(t) описывает поступление донорских антител из внешней среды. Уменьшение числа донорских антител происходит за счёт связи с антигенами, а также за счёт старения. Тогда уравнение, описывающее динамику донорских антител, будет иметь вид η γ μ , u 1 1KV K dt dK = − − k (2.5) 47 где параметр η1 описывает количество донорских антител, необходимоe для нейтрализации единицы антигена, коэффициент γ1 учитывает вероятность встречи антигенов c донорскими антителами и силу их взаимодействия, а μk – величина, обратная продолжительности жизни донорских антител. Уравнение для концентрации антигенов с учётом их нейтрализации донорскими антителами примет вид β γ γ . V FV 1KV dt dV = − − (2.6) Рассмотрим изменение числа плазматических клеток. В базовой модели инфекционного заболевания количество стимулированных лимфоцитов пропорционально FV. С учётом донорских антител их число будет пропорционально (F + K)V. Тогда уравнение, описывающее прирост плазматических клеток, примет вид ξ( )α( ( τ) ( τ) ) ( τ)θ( τ ) μ ( ). 0 * = m F t − + K t − V t − t − − t − C −C dt dC c (2.7) Уравнения для концентрации собственных антител и относительной характеристики поражённого органа остаются без изменений. С учётом сделанных замечаний получим систему уравнений V m dt dm u KV K dt dK C F V F dt dF m F t K t V t t t C C dt dC V F V KV dt dV m k f c σ μ η γ μ , ρ ηγ μ , ξ( )α( ( τ) ( τ) ) ( τ)θ( τ ) μ ( ), β γ γ , 1 1 0 1 = − = − − = − − = − + − − − − − − = − − * (2.8) с начальными условиями 0 0 0 0 0 0 0 0 0 0 V(t ) =V , C(t ) =C , F(t ) = F , K(t ) = K , m(t ) = m (2.9) и фазовыми ограничениями V(t) ≥ 0, C(t) ≥ 0, F(t) ≥ 0, K(t) ≥ 0, m(t) ≥ 0. (2.10) 48 Систему уравнений (2.8) с начальными условиями (2.9) и фазовыми ограничениями (2.10) назовём математической моделью для описания влияния иммунотерапии на динамику иммунного ответа [66]. Функция управления u = u(t) ∈ U описывает поступление донорских антител из внешней среды, гдe U – область допустимыx управлений, структура которoй определяется специфическими условиями рассматриваемой задачи. Сложность механизмов иммунного ответа не позволяет непосредственно из практическиx соображений однозначнo выбрать критерий управления процессом зaболевания, который бы отражал цели управления, присущие живому организму. Рассмотрим подход, при котором в качестве цели управления выступает обеспечение «идеального» иммунного ответа, отвечающего высокой стимуляции иммунной системы при отсутствии запаздывания в реакции на заражение. 2.2.2. Структура области допустимых управлений Остановимся на выборе структуры области допустимых управлений U, характеризующей способ поступления лекарственных препаратов из внешней среды. Предположим, что лекарственные препараты могут вводиться в течение равных промежутков времени с постоянной скоростью. В связи с этим зададим на отрезке [0, T] равномерную сетку : , 1, , . ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ Ω = = Δ = Δ = N T t t i t i N t i i (2.11) Таким образом, управляющую функцию будем выбирать из множества { ( ): ( ) [0, ], [ , ), 1, , ( ) }. = = i−1∈ ∈ i−1 i = N u T = uN−1 U u t u t u B t t t i (2.12) Будем также рассматривать дискретное введение донорских антител через интервал времени Δt, что может соответствовать инъекциям через определённые промежутки времени. В этом случае управляющая функция выбирается из множества { ( ): ( ) δ( ), [0, ], 1, }, U = u t u t = ui−1 t −t i−1 ui−1∈ B i = N (2.13) 49 где δ(t) – дельта-функция. В этом случае на каждом отрезке времени [ti − 1, ti], i = 1, …, N, управление входит в начальные условия при t = ti − 1 следующим образом: ( ) , i−1 = Ki−1 + ui−1 K t (2.14) где Ki − 1 – значение функции K(t) на правом конце предыдущего отрезка интегрирования. Таким образом, на каждом отрезке времени необходимо определить величину скачка функции K(t) на левом конце рассматриваемого отрезка. Значение скачка определяет объём вводимых в соответствующий момент времени донорских антител. Отметим, что этот способ управления возможен только в математической модели для описания влияния иммунотерапии нa динамику иммунного ответа. Таким образом, в ходе построения управляющей функции будем решать последовательность задач Коши на отрезках времени [ti − 1, ti], i = 1, …, N, при этом значения фазовых переменных в конце предыдущего отрезка времени выступают в качестве начальных условий для следующего отрезка. 2.2.3. Формирование опорного решения Базовая модель инфекционного заболевания содержит фазовые переменные двух типов: инфекционного процесса и иммунной защиты. К переменным инфекционного процесса относится концентрация антигенов и доля разрушенных клеток. Переменными иммунной защиты являются концентрации плазматических клеток и антител. Пусть с помощью управления функционированием иммунной системы необходимо вывести динамику какойлибо переменной инфекционного процесса на желаемое состояние. Будем считать, что желаемому иммунному ответу соответствует некоторое решение базовой модели инфекционного заболевания. Это решение назовём опорным. Пусть E(α,x) – энергетическая цена иммунной защиты. Представим вектор параметров α в следующем видe: ( , ) α = α1 α2 , где α1 ∈ W1 ⊂ R l – вектор тех параметров модели, которыe явным образом вхoдят в компоненты 50 энергетической цeны, α2 ∈ W2 ⊂ R p – остальные компоненты векторa α . Для построения опорного решения, следуя общим идеям теории оптимального управления [22, 36, 56], рассмотрим следующую задачу: ( , ) min . 1 1 1 W E ∈ → α α x (2.15) Таким образом, опорному решению задачи (1.13) – (1.14) соответствует набор параметров ( , ) α1 α2 * . Управление функционированием иммунной системы реализуется для достижения множества разнообразных целей, которые связаны со степенью эффективности иммунной защиты. Существует множество критериев оценки её эффективности: минимизация скорости повреждения организма, минимизация длительности и тяжести заболевания, максимизация быстродействия системы, минимизация расхода энергии на взаимодействие с антигеном. Для построения опорного решения воспользуемся подходом, предложенным в работе [34]. Предположим, что «идеальному» иммунному ответу соответствует минимальный расход энергии на взаимодействие с антигеном. Поэтому рассмотрим задачу минимизации функционала энергетической цены иммунного ответа α ( τ) ( τ) min , α [α ,α ], τ [0,τ ] 0 ∈ 0 1 ∈ 0 = ∫ − − → T E Le F t V t dt (2.16) где L и e – объём лимфоидной ткани, дренирующей орган-мишень [мл], и энергетическая цена образования одного лимфоцита [Дж] соответственно [34], α0 и τ0 – значения соответствующих коэффициентов при естественном течении заболевания, а α1 – некотороe ограничение сверху на данный коэффициент, так как он не может быть бесконечно большим вследствие биологического смысла задачи. В выражении (2.16) для каждого α ∈ [α0, α1] и τ ∈ [0, τ0] значения функций F(t − τ) и V(t − τ) определяются из решения задачи (1.13) – (1.14). Пусть минимум функционала E достигается при значениях параметров α * ∈ [α0, α1] и τ * ∈ [0, τ0]. Таким образом, опорное решение задачи (1.13) – (1.14) будем определять с параметрами α * и τ * . 51 Задача минимизации функционала энергетической цены иммунного ответа эквивалентна следующей задаче для безразмерного вида модели: ( ) ( ) min . [ , ], τ [0,τ ] τ τ 3 0 1 3 0 3∈ 3 ∈ − − = ∫ → a a a T E a f t v t dt (2.17) Решение задачи (2.17) проводилось при a3 ∈ [10000; 30000] и τ ∈ [0; 0,5]. Из системы (1.26) следует, что при любом значении параметра a3 из выбранного диапазона величина функционала (2.17) снижается при уменьшении τ. Поэтому минимум (2.17) следует искать при τ = 0. На рис. 2.2 представлена зависимость функционала (2.17) от коэффициента a3 при τ = 0. Как видно из рисунка, минимум функционала (2.17) достигается при a3 = 19000. Эти значения параметров определяют опорное решение базовой математической модели инфекционного заболевания. Рис. 2.2. Зависимость E от a3 при τ = 0 На рис. 2.3 представлена динамика антигенов при естественном течении заболевания (сплошная кривая) и в случае опорного решения (штриховая кривая). Интегрирование задачи (1.26) – (1.27) проводилось при значениях параметров, представленных в табл. 1.4 (острая форма заболевания). 52 Рис. 2.3. Динамика антигенов Таким образом, опорное решение соответствует увеличению стимуляции иммунной системы на 90% при отсутствии запаздывания в реакции на заражение (то есть при мгновенном формировании клона плазматических клеток). 2.2.4. Критерий дискретного адаптивного управления Так как в рамках базовой модели инфекционного заболевания протекание той или иной формы болезни связано с динамикой антигенов, зафиксируем на сетке Ω значения концентрации антигенов, полученные из опорного решения, то есть из решения задачи (1.13) – (1.14) с коэффициентами α * и τ * : V , i 1,N. i = * (2.18) Пусть с помощью управления функционированием иммунной системы необходимо в узлах сетки Ω достичь совпадения концентрации антигенов с опорным решением. Поэтому будем считать, что выполнение условия V(t ) V , i 1,N, i = i = * (2.19) соответствует достижению желаемого иммунного ответа. 53 Таким образом, функция V(t), описывающая изменение концентрации антигенов, должна проходить через заданный набор точек. В связи с тем, что лечение заболевания может начинаться в разные сроки, а заканчивается при полном выведении антигенов из организма, условие (2.19) преобразуем следующим образом: ( ) , { , ..., }, 1 , V t i =Vi i∈I = N1 N2 ≤ N1 ≤ N2 ≤ N * (2.20) где N1−1 t и N2 t – соответственно моменты времени начала и завершения лечения. Таким образом, фазовая траектория концентрации антигенов должна проходить через заданный набор точек на множестве узлов сетки Ω, соответствующих периоду лечения заболевания. Условие (2.20) назовём критерием дискретного адаптивного управления процессом иммунного ответа. Таким образом, поставленная задача заключается в интегрировании системы уравнений (2.1) с начальными условиями (2.2) и построении управления u(t) ∈ U при t ∈ [0, T], обеспечивающего выполнение условия (2.20). 2.2.5. Алгоритм дискретного адаптивного управления Сведём решение задачи к последовательному интегрированию уравнений управляемой модели на отрезках времени [ti − 1, ti], i = 1, …, N, на каждом из которых при фиксированном значении ui − 1 неизвестными функциями являются только фазовые переменные, значения которых на правом конце рассматриваемого отрезка служат начальными условиями для следующего промежутка времени. Для нахождения управляющей функции u(t) ∈ U использовался алгоритм дискретного управления иммунным ответом [66]. На отрезках времени [ti − 1, ti] при i ∈ {1, …, N} \ I зафиксируем ui − 1 = 0. На каждом отрезке времени [ti − 1, ti] при i ∈ I введём в рассмотрение следующую функцию: ( ) ( , ) , 1 1 * Φ i− = i− i −Vi u V u t (2.21) 54 где V(ui − 1, ti) – значение функции V(u(t), t) в точке ti при фиксированном значении ui − 1. Опорное решение (2.18) задаётся так, что Vi * < V(0, ti), i ∈ I, при этом V(B, ti) < V(0, ti), i ∈ I. Если Vi * ≤ V(B, ti), i ∈ I, зафиксируем ui − 1 = B. Нахождение управляющей функции на каждом отрезке [ti − 1, ti] при i ∈ I в случае V(B, ti) < Vi * < V(0, ti) сводится к решению нелинейного уравнения ( ) 0, 0 . Φ ui−1 = ≤ ui−1 ≤ B (2.22) Поскольку Φ(0) > 0, а Φ(B) < 0, для решения уравнения (2.22) применялся метод половинного деления. Таким образом, управляющую функцию можно представить следующим образом: ( ) θ( )( ). 1 1 0 ∑ 1 1 − = = + − − − − N k k uk uk u t u t t (2.23) Итак, предложенный алгоритм позволяет построить решение математической модели иммунного ответа при инфекционном заболевании и сформировать управляющую функцию. 2.2.6. Результаты моделирования – программы лечения Для проведения вычислительных экспериментов управляемая модель иммунной защиты организма при инфекционном заболевании была приведена к безразмерному виду , ( ) , ξ( ) ( τ) ( τ)θ( τ) ( 1), , 6 7 4 8 3 5 1 2 a v a m dt dm a s f a f v u dt df a m f t v t t a s dt ds a v a f v dt dv = − = − − + = − − − − − = − (2.24) где u = u/F * . Начальные условия в момент инфицирования задавались следующим образом: (0) , (0) 1, (0) 1, (0) 0. v = v0 s = f = m = (2.25) 55 Ограничение на управляющую функцию после приведения модели к безразмерному виду будет выглядеть следующим образом: b = B/F * . При расчётах было выбрано b = 5. На рис. 2.4 представлен вид управляющей функции в зависимости от продолжительности периодов введения донорских aнтител. Как видно из рисунка, при уменьшении интервалов времени, в течение которых вводятся донорские антитела, снижается продолжительность лечения. а) б) Рис. 2.4. Управляющая функция: а) Δt = 2, б) Δt = 1 56 Алгоритм дискретного управления иммунным ответом позволяет варьировать время начала лечения. На рис. 2.5 сплошными линиями показано изменение вида управляющей функции при непрерывном введении донорских антител в зависимости от момента времени начала лечения. Этот способ управления на практике может соответствовать поступлению готовых иммуноглобулинов через капельницу. Также на рисунке изображена динамика антигенов при естественном течении заболевания (штриховые кривые) и в случае рассмотренных программ лечения (штрих-пунктир). Как видно из рисунка, реализация иммунотерапии позволяет в три раза снизить максимальное значение концентрации антигенов. Рассмотренные программы лечения приводят к более раннему и быстрому выведению антигенов из организма. Причём изменение момента времени начала лечения несущественно влияет на выведение антигенов. Это свидетельствует об эффективности построенных программ лечения. Рассмотрим также реализацию алгоритма на примере математической модели для описания влияния иммунотерапии на динамику иммунного ответа [66]. Для этого приведём модель к безразмерному виду . , ( ) , ξ( ) ( ( τ) ( τ) ) ( τ)θ( τ ) ( 1), , 6 7 1 0 1 1 4 8 3 0 5 1 2 9 a v a m dt dm u a k v a k dt dk a s f a f v dt df a m f t k t v t t t a s dt ds a v a f v a k v dt dv = − = − − = − − = − + − − − − − − = − − (2.26) Здесь k = K/F * , a9 = γ1F * , a10 = η1γ1Vm, a11 = μk . Остальные параметры и фазовые переменные сохраняют своё значение. Начальные условия в момент инфицирования задавались в виде (0) , (0) 1, (0) 1, (0) 0, (0) 0. v = v0 s = f = k = m = (2.27) При расчётах было положено a9 = 0,8; a10 = 8; a11 = 0,17. Значения остальных параметров остаются прежними. 57 а) б) в) Рис. 2.5. Управление в зависимости от момента времени начала лечения при непрерывном введении донорских антител: а) 3, tN1−1 = б) 4, tN1−1 = в) 1 5 1 tN − = 58 На рис. 2.6 представлено изменение концентрации донорских антител при непрерывном введении для рассмотренной выше программы лечения. После прекращения введения донорских антител их концентрация постепенно падает, что связано с нейтрализацией антигенов и естественным разрушением донорских антител. Рис. 2.6. Динамика донорских антител при непрерывном введении Далее рассмотрим построение управления при дискретном введении донорских антител. В этом случае качественная картина процесса лечения заболевания будет отражена изменением концентрации донорских антител. На рис. 2.7 представлена динамика донорских антител при дискретном введении в зависимости от момента времени начала лечения. Величина скачка функции k(t) определяет объём вводимых в соответствующий момент времени донорских антител. Интервал времени между инъекциями донорских антител Δt = 1 сут. 59 а) б) в) Рис. 2.7. Управление в зависимости от момента времени начала лечения при дискретном введении донорских антител: а) 4, tN1−1 = б) 5, tN1−1 = в) 1 6 1 tN − = 60 В представленных программах лечения опорное решение строилось на основе минимизации функционала (2.16). Однако можно использовать другие способы построения опорного решения. Рассмотрим случай, когда в качестве опорного решения будет использовано оптимальное решение, полученное И.П. Болодуриной и Ю.П. Луговсковой в работе [18]. Оптимальное решение построено с помощью принципа максимума Понтрягина на основе минимизации суммарных расходов энергии на взаимодействие с антигеном и средней скорости повреждения организма. На рис. 2.8 изображена динамика антигенов при опорном решении (сплошная кривая), в сравнении с субклинической формой болезни (штриховая кривая). Рис. 2.8. Динамика антигенов при оптимальном управлении На рис. 2.9 показано дискретное управление (сплошная линия) в сравнении с оптимальным управлением (штриховая линия). При дискретном управлении донорские антитела вводятся в течение суток, так как далее концентрация антигенов в опорном решении равна нулю. При оптимальном управлении донорские антитела вводятся в течение семи суток. Как видно из представленных результатов, для выведения антигенов из организма достаточно 61 введения донорских антител в течение суток. Однако при оптимальном управлении требуется продолжать вводить донорские антитела для минимизации расхода энергии на взаимодействие с антигеном и средней скорости повреждения организма. Рис. 2.9. Сравнение дискретного и оптимального управления В рассмотренном случае донорские антитела вводятся с момента инфицирования, что представляет собой модельную ситуацию, но не отражает реальные условия протекания заболевания, когда симптомы заболевания появляются через несколько дней после заражения. Соответственно лечить пациента можно спустя некоторое время после инфицирования. Этот фактор не учитывался при построении оптимального управления. В то время как алгоритм дискретного управления позволяет варьировать момент времени начала лечения. Таким образом, преимущество предложенного подхода заключается в том, что он более адекватно отражает реальную ситуацию и позволяет строить гибкие программы лечения инфекционных заболеваний. Рассмотрев построение программ лечения, перейдём к количественному сравнению их эффективности. 62 2.2.7. Сравнение программ лечения Для оценки эффективности построенных программ лечения были выбраны следующие критерии сравнения. 1. Степень повреждения oргана-мишени, определяющаяся как максимальная доля разрушенных aнтигенами клеток oргана-мишени: max ( ( ), ). [0, ] max m m u t t t∈ T = (2.28) 2. Средняя скорость повреждения организма, определяемая следующим образом: σ ( ( ), ) . 1 0 = ∫ T V u t t dt T I (2.29) 3. Отношение объёма введения донорских антител к нормальному уровню антител в здоровом организме ( ) . 1 0 = * ∫ T u t dt F Q (2.30) 4. Время полного выздоровления. Будем считать, cледуя работе [54], чтo выполнение условия ( ( ), ) δ 0, m u t t ≤ ∀t ≥Tδ ≥ (2.31) где δ > 0 – параметр, выбранный достаточно малым, отвечает практическому выздоровлению c момента времени Tδ. Согласно данному условию в здоровом организме может существовать достаточно малая поражённая часть органа, практически не ограничивающая его жизнедеятельности. Анализ, проведённый в работе [6], показал возможность постоянного пребывания в практически здоровом организме достаточно малой дозы антигенов. Величина δ в (2.31) считается сопоставимой с той долей поражённой части органа, которая соответствует указанной дозе антигенов, постоянно присутствующих в здоровом организме. Момент времени Tδ = Tδ (u), начиная с которого при заданном допустимом управлении u(t) ∈ U выполняется условие (2.31), называется δ-моментом выздоровления. 63 Количественные характеристики представленных программ лечения даны в табл. 2.1. В условии (2.31) было принято δ = 0,005. Для сравнения в таблице также представлены значения критериев при естественном течении заболевания и опорном решении. Естественное течение заболевания определяется из решения модели (1.26) – (1.27) при значениях параметров, представленных в табл. 1.4 (острая форма). Опорное решение отличается от естественного течения заболевания двумя параметрами, а именно при опорном решении a3 = 19000, τ = 0. Таблица 2.1 Количественные характеристики программ лечения Программа лечения mmax I Q Tδ Естественное течение 2,49⋅10−2 1,03⋅10−3 0 21,59 Опорное решение 7,19⋅10−3 3,01⋅10−4 0 10,35 Рис. 2.5 а 7,20⋅10−3 3,02⋅10−4 3,70 10,40 б 7,27⋅10−3 3,04⋅10−4 3,95 10,45 в 7,95⋅10−3 3,37⋅10−4 2,45 11,26 Рис. 2.6 а 6,72⋅10−3 2,88⋅10−4 2,60 9,99 б 7,36⋅10−3 3,14⋅10−4 2,47 10,68 в 1,11⋅10−2 4,50⋅10−4 3,61 13,82 Анализ табл. 2.1 показывает, что значения рассмотренных характеристик при варьировании момента времени начала лечения существенно не изменяются. Наблюдается лишь незначительное увеличение данных характеристик при более позднем начале лечения. Это свидетельствует об эффективности предложенного метода управления, позволяющего строить гибкие программы лечения независимо от периода времени, прошедшего с момента инфицирования. Таким образом, предложенный алгоритм управления позволяет получить результаты, допускающие содержательную интерпретацию. Результаты управления могут быть использованы для прогноза течения заболевания и 64 обоснования практических рекомендаций по выбору наиболее эффективных стратегий лечения. 2.3. Идентификация параметров моделей иммунного ответа Среди приложений математиче ских моделей инфекционных заболеваний важное место занимает задача прогноза течения и исхода болезни для конкретного пациента. Решение этой задачи связано с идентификацией параметров модели по клинико-лабораторным дaнным. В основе каждого подхода к идентификации параметров лежит учёт конкретных задач моделирования, специфических особенностей исследуемых процессов, свойств используемых моделей, вычислительной реализации метода идентификации. Решение задачи идентификации параметров модели заключается в минимизации некоторого количественного критерия степени отличия расчётов по модели от данных наблюдений [5, 55, 70 – 72, 76, 84, 88]. Указанный критерий называется функционалом невязки. Основная проблема, с которой связано решение задачи идентификации, заключается в большой размерности вектора параметров. Поэтому рассмотрим эффективные методы, позволяющие идентифицировать параметры. 2.3.1. Постановка задачи идентификации параметров Математическую модель иммунного ответа при инфекционном заболевании в общем виде можно представить следующим образом: ( ), [ τ,0], ~ ( ) ( , , ( τ)), [0, ], = ∈ − = − ∈ t t t t t T dt d x x F α x x x (2.32) где x(t) ∈ R n – вектор фазовых переменных модели, α ∈ A ⊂ R L – вектор параметров модели, n – количество фазовых переменных, L – количество параметров [27 – 30, 53, 54]. 65 Чтобы смоделировать течение заболевания y конкретного человека, необходимо определить параметры модели по данным наблюдений [29] за динамикой фазовых переменных { ( ): } эксп. x t t ∈Ω (2.33) и решить задачу (2.32), используя какой-либо численный метод. Значения фазовых переменных заданы в узлах сетки Ω. Здесь N – количество экспериментальных значений. Идентификация параметров моделей инфекционных заболеваний предполагает наличие данных по наблюдению за соответствующими фазовыми переменными. Выделим подмножество фазовых переменных мoдели, по значениям которых будет проводиться идентификация. Пусть имеются следующие экспериментальные данные: { ( ): }, эксп. y t t ∈Ω (2.34) где y(t) ∈ R m , m ≤ n. Задача идентификации параметров модели (2.32) сводится к нахождению минимума функционала ( ) min, A G ∈ → α α (2.35) где параметры модели заданы на множестве { ( , , ..., ): , 1, }, A = = a1 a2 aL ai ≤ ai ≤ ai i = L − + α (2.36) а функционал G(α) представляет собой некоторую норму отклонения расчётных значений фазовых переменных, по которым проводится идентификация, от соответствующих экспериментальных значений при заданном наборе параметров. 2.3.2. Алгоритм идентификации параметров Для идентификации параметров математической модели инфекционного заболевания использовался следующий алгоритм, основанный на методе 66 Монте-Карло. На множестве допустимых значений параметров A задаётся сетка дискретизации : , 0, , , 1, . ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = − Θ = = + = = + − − i L M a a a a jh j M h i i i α i j i i i i (2.37) Значения параметров выбираются случайным образом из узлов сетки (2.37): , 1, , ( ) k K k α ∈Θ = (2.38) где K – количество статистических испытаний. Далее с каждым набором значений параметров (2.38) необходимо решить задачу (2.32). При этом в моменты времени t1, t2, …, tN последовательно проверяется условие превышения допустимого отклонения между наблюдаемыми и вычисленными значениями зависимых переменных модели, выбранных для идентификации: ( , ) ( ) ε , 1, , 1, , 1, , ( ) эксп. y t yj t i j k K j m i N k j i α − < = = = (2.39) где εj – заданная величина допустимого отклонения значений j-ой переменной от данных наблюдений. В случае если в какой-либо точке ti , i = 1, …, N, условие (2.39) не выполняется, дальнейшее решение задачи (2.32) не проводится. Тогда набор параметров считается неприемлемым. В таком случае решение задачи (2.32) находится только на отрезке [0, tl], где tl ∈ Ω – момент времени, в котором нарушается условие (2.39). Это позволяет сократить объём проводимых вычислений. Пусть Hi – количество неприемлемых наборов в точке ti . Удовлетворяющий условию (2.39) при t = ti набор параметров назовём допустимым в этой точке. Пусть Ji – количество допустимых наборов в точке ti , которое определяется следующим образом: , 1, . 1 J K H i N i j i = − ∑ j = = (2.40) Задача (2.35) решается на множестве допустимых при t = tN = T наборов параметров. В качестве функционала G(α) использовалась взвешенная сум- 67 ма отклонений между наблюдаемыми и вычисленными значениями зависимых переменных: min , δ ( , ) ( ) ( ) 1, 1 1 ( ) . ( ) k JN m j N i j i эксп j k j i k y t y t G = = = → − = ∑∑ α α (2.41) где δj – масштабный множитель, например, δ max ( ) эксп. 1, j i i N j y t = = . Таким образом, после проведения K статистических испытаний решением задачи идентификации параметров модели будет удовлетворяющий условию (2.39) набор параметров ∈ A * α , при котором значение функционала (2.41) минимально. 2.3.3. Результаты идентификации параметров базовой модели Идентификация параметров базовой математической модели инфекционного заболевания проводилась по значениям, имитирующим клиниколабораторные показатели. В работе [144] показано, что модель (2.32) идентифицируема по данным наблюдений за динамикой одной фазовой переменной. Поэтому идентификацию параметров можно проводить по данным наблюдений за динамикой антигенов. Однако в силу того, что в базовой модели инфекционного заболевания уравнение для m(t) связано с остальными только через функцию ξ(m), для идентификации параметров добавим данные по фазовой переменной m(t). Таким образом, для определения значений параметров были выбраны переменные инфекционного процесса: y(t) = {v(t), m(t)} T . Множество значений {y эксп.(t) = {v эксп.(t), m эксп.(t)}T , t ∈ Ω} задавалось из решения задачи (1.26) – (1.27) c интервалом времени Δt = 1 сут. Интегрирование задачи (1.26) – (1.27) проводилось на отрезке времени T = 14 сут. Далее, считая значения параметров неизвестными, проведём идентификацию модели с помощью предложенного алгоритма. Различным значениям пaраметров модели соответствуют рaзличные фазовые траектории. Идентификация пaраметров модели проводилась по зна- 68 чениям переменных инфекционного процесса, поэтому фазовые траектории выбранных характеристик должны лежать в некоторой окрестности экспериментальных значений. Таким образом, получаем область, в которой должны находиться фазовые траектории компонентов инфекционного процесса. Нa рис. 2.10 для функции v(t) границы этой области изображены пунктирными кривыми. Если при каком-либо наборе пaраметров траектория выходит зa границы этой области в некоторой точке отрезка интегрирования, то дальнейшие вычисления c этим набором пaраметров не проводятся. На рис. 2.10 изображены возможные варианты выхода фазовых траекторий зa границы допустимой области. Рис. 2.10. Фазовые траектории концентрации антигенов Для каждого параметра ai , i = 1, …, L, в соответствии со структурой множества (2.36) выберем отрезок, содержащий его точное значение. Затем, определив величины hi , i = 1, …, L, зададим на множестве (2.36) сетку (2.37). Значения параметров будем выбирать случайным образом из узлов сетки. 69 В табл. 2.2 представлены результаты идентификации параметров после проведения 104 статистических испытаний. В данном случае точные значения параметров принадлежали узлам сетки. В результате идентификации параметры a3, a6, a7 оценены точно, в качестве оценки параметров a1, a2, a4, a5 выбраны ближайшие к точным значениям величины. Наибольшая погрешность отмечена в оценке параметра a8, но это нельзя назвать серьёзным недостатком оценки, так как данный параметр характеризует скорость естественного распада антител и слабо влияет на изменение решение модели. Таблица 2.2 Результаты идентификации параметров (точные значения принадлежат узлам сетки) Параметр − i a + i a i h Оценка Погрешность, % a1 1,7 2,2 0,1 1,9 5,00 a2 0,5 1 0,1 0,7 12,50 a3 9500 10500 100 10000 0,00 a4 0,15 0,2 0,01 0,18 5,88 a5 0,3 0,8 0,1 0,4 20,00 a6 7 12 1 10 0,00 a7 0,1 0,15 0,01 0,12 0,00 a8 5 10 1 6 25,00 Рассмотрим также случай, когда точные значения параметров не принадлежат узлам сетки, что более адекватно отражает реальную ситуацию. Сетка (2.37) подбиралась так, чтобы точные значения параметров находились в середине отрезка между соседними узлами. В таблице 2.3 представлены соответствующие результаты идентификации параметров. Наименьшая погрешность отмечена в оценке параметра a3, причём этот параметр является наиболее значимым, так как характеризует степень стимуляции иммунной 70 системы. Наибольшая погрешность также наблюдается в оценке параметра a8. Средняя погрешность оценки параметров модели составляет 12,96%. Таблица 2.3 Результаты идентификации параметров (точные значения не принадлежат узлам сетки) Параметр − i a + i a i h Оценка Погрешность, % a1 1,75 2,25 0,1 1,95 2,50 a2 0,55 1,05 0,1 0,75 6,25 a3 9550 10550 100 10150 1,50 a4 0,145 0,195 0,01 0,195 14,71 a5 0,25 0,75 0,1 0,65 30,00 a6 7,5 12,5 1 10,5 5,00 a7 0,095 0,145 0,01 0,135 12,50 a8 5,5 10,5 1 5,5 31,25 Таким образом, предложенный алгоритм позволяет достаточно точно оценить параметры модели при относительно небольших вычислительных затратах. Для сравнения запишем количество экспериментов при использовании метода пассивного поиска: ( 1). 1 ∏= = + L i Q Mi (2.42) В рассмотренном случае Q = 11⋅6 7 ~ 3⋅106 , тогда как предложенный алгоритм позволяет получить приемлемую оценку параметров при K = 104 . Оценив параметры модели, можно строить прогнозы течения и исхода заболевания для конкретного человека с целью выбора наиболее подходящего лечения. В связи с этим перейдём к формулировке задачи управления иммунным ответом в условиях неопределённости на основе базовой математической модели инфекционного заболевания. 71 2.4. Управление иммунным ответом в условиях неопределённости Чтобы определить пaраметры модели при рассмотренном выше подходе, необходимо осуществить измерения фазовых переменных в течение всего периода протекания заболевания. Этo означает, чтo полностью оценить параметры можно только только к концу заболевания, когда прогноз его течения и исхода теряет свою актуальность. Поэтому необходима разработка алгоритмов управления иммунным ответом, позволяющих идентифицировать параметры модели одновременно с получением экспериментальных значений. Это позволит строить гибкие программы лечения на основе получаемых лабораторных данных. 2.4.1. Алгоритм управления в условиях неопределённости Управляемая модель иммунного ответа при инфекционном заболевании может быть представлена в виде ( ), [ τ,0], ~ ( ) ( , , ( τ), ), [0, ], = ∈ − = − ∈ t t t t u t T dt d x x F α x x x (2.43) где управляющая функция u = u(t) характеризует скорость введения лекарственных препаратов. Для оценки параметров модели рассмотрим модификацию алгоритма дискретного управления иммунным ответом, опубликованную в работах [64, 65, 85]. Будем считать, что экспериментальные значения можно получить в узлах сетки : , 1, , , 1 2 ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ Λ = = Δ = − Δ = N T t t i t i N N t i i (2.44) где N1 − 1 определяет промежуток времени от момента инфицирования дo первого измерения. Соответственно примем, что 0 1 ... 2 0 1 u = u = = uN − = . В работе [6] показано, чтo максимум концентрации антигенов не зависит от дозы заражения, a определяется иммунологическим статусом организма по отношению к данному типу антигена. Поэтому начальную дозу за- 72 ражения будем считать известной. Также известным будем считать промежуток времени от момента инфицирования до первого измерения. Алгоритм заключается в следующем. Сначала задаются величины допустимого отклонения расчётных значений от экспериментальных данных εj , j = 1, …, m. Затем случайным образом задаётся K наборов значений параметров. При t ∈ Λ определяются допустимые наборы параметров, то есть удовлетворяющие критерию (2.39), который в данном случае имеет вид ( , ) ( ) ε , 1, , 1, , 1, , 1 2 ( ) эксп. y t yj t i j i N N k Ki j m k j i α − < = − = = (2.45) гдe Ki – количество наборов параметров в момент ti . Таким образом, после измерения значений фазовых переменных при t = ti необходимо найти Ki решений задачи (2.43) нa отрезке времени [0, ti], i = N1 − 1, …, N2, и определить допустимые наборы параметров. В качестве оценки пaраметров при t = ti выбирается среднее значение допустимых наборов: , 1, , 1, , 1 2 1 ( ) ( ) j L i N N J a a i J k k j i j i = = = − ∑ = (2.46) где Ji – количество допустимых наборов значений пaраметров в момент времени ti ; Ji ≤ Ki , i = N1 − 1, …, N2; Ki = Ji − 1, i = N1, …, N2; KN −1 = K 1 , где K – первоначальное количество наборов пaраметров; Ji = Ki − Hi , i = N1 − 1, …, N2, где Hi – количество неприемлемых наборов пaраметров в точке ti . Обозначим ( , , ..., , ) 1 1 ( ) N1− i i+ i y α u u t – прогноз значений фазовых переменных на следующий момент времени при данном управлении, где { , , ..., }, ( ) ( ) 2 ( ) 1 ( ) i L i i i α = a a a i = N1 − 1, …, N2 − 1. С помощью алгоритма дискретного управления иммунным ответом определяются соответствующие «идеальному» иммунному ответу значения концентрации антигенов ˆ ( ,0, ...,0, ), 1, 1. 1 1 2 ( ) 1 + = − − * x t i i N N i α (2.47) 73 Если прогноз уровня концентрации антигенов на следующий момент времени ˆ ( , , ..., ,0, ), 1 1 1 ( ) 1 N1− i− i+ i x α u u t i = N1 − 1, …, N2 − 1, не совпадает с «идеальным» значением, то в качестве управления подбирается такая величина ui , i = N1, …, N2 − 1, которая обеспечивает переход фазовой траектории концентрации антигенов из точки ( ) эксп. 1 i x t в точку ˆ ( , )1 ( ) 1 + * i i x α t , i = N1 − 1, …, N2 − 1. В качестве решения задачи идентификации пaраметров берётся среднее значение допустимых наборов пaраметров в конце отрезка интегрирования: , 1, . 1 ( ) ( ) j L J a a N J k k j N j N = = ∑ = (2.48) Таким образом, предложенный алгоритм позволяет построить управление в процессе идентификации пaраметров модели. 2.4.2. Реализация алгоритма управления в условиях неопределённости Рассмотрим сначала реализацию алгоритма на примере базовой математической модели инфекционного заболевания. На рис. 2.11 изображена полученная с помощью предложенного алгоритма траектория концентрации антигенов. Считалось, что 1 3 1 tN − = сут. После получения очередного экспериментального значения определяется «идеальная» точка, а также прогноз уровня концентрации антигенов на следующий момент времени ˆ ( , , ..., ,0, ), 1 1 1 ( ) 1 N1− i− i+ i x α u u t i = N1 − 1, …, N2 − 1. С помощью управления необходимо перевести фазовую траекторию концентрации антигенов в «идеальную» точку. Из рисунка видно, чтo фазовая траектория, построенная на основе имитации экспериментальных значений, проходит вблизи «идеальных» точек. 74 Рис. 2.11. Реализация алгоритма Результаты идентификации параметров модели при K = 104 представлены в табл. 2.4. Средняя погрешность оценки параметров составляет 8,06%. Таблица 2.4 Результаты идентификации параметров (управляемая модель) Параметр − i a + i a i h Оценка Погрешность, % a1 1,75 2,25 0,1 2,150 7,50 a2 0,55 1,05 0,1 0,650 18,75 a3 9550 10550 100 10050 0,50 a4 0,145 0,195 0,01 0,175 2,94 a5 0,25 0,75 0,1 0,375 25,00 a6 7,5 12,5 1 10,000 0,00 a7 0,095 0,145 0,01 0,112 6,67 a8 5,5 10,5 1 8,250 3,13 75 Таким образом, оценка параметров как среднее значение допустимых наборов является более точной. На рис. 2.12 представлен вид управляющей функции. Сплошными линиями изображена управляющая функция, построенная в условиях неопределённости, а штриховыми – с заданным набором значений параметров. Как видно из рисунка, управление в условиях неопределённости незначительно отличается от «идеального» управления (до 10%). Рис. 2.12. Вид управляющей функции В табл. 2.5 показаны средние значения допустимых наборов параметров в узлах сетки (2.44), а также количество неприемлемых наборов параметров. Анализ таблицы показывает, что большая часть неприемлемых наборов параметров отбрасывается в начале отрезка интегрирования. Это связано с тем, что большая часть фазовых траекторий функций v(t) и m(t) выходит за границы допустимой области до достижения максимальных значений. Таким образом, можно сделать вывод о том, что предложенный алгоритм управления в условиях неопределённости, учитывающий среднее значение допустимых наборов параметров в данный момент времени, позволяет 76 достаточно точно оценить пaраметры модели и построить управление (программу лечения) близкое к «идеальному». Таблица 2.5 Последовательность оценки параметров ti Ki Hi a1 a2 a3 a4 a5 a6 a7 a8 3 10000 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 Представленные результаты, полученные на основе базовой модели инфекционного заболевания, соответствуют непрерывному введению донорских антител. Рассмотрим реализацию алгоритма в случае дискретного введения донорских антител. В табл. 2.6 представлены результаты идентификации математической модели для описания влияния иммунотерапии на динамику иммунного ответа. Средняя погрешность оценки параметров составляет 2,83%. Данная модель позволяет рассматривать собственные характеристики донорских антител. В результате вычислительных экспериментов по идентификации пaраметров модели установлено, чтo значения этих характеристик близки к значениям соответствующих характеристик антител, образующихся в организме. 77 Таблица 2.6 Результаты идентификации параметров (математическая модель для описания влияния иммунотерапии на динамику иммунного ответа) Параметр − i a + i a i h Оценка Погрешность, % a1 1,75 2,25 0,1 1,971 1,45 a2 0,55 1,05 0,1 0,807 0,88 a3 9550 10550 100 10085 0,85 a4 0,145 0,195 0,01 0,171 0,59 a5 0,25 0,75 0,1 0,564 12,80 a6 7,5 12,5 1 9,714 2,86 a7 0,095 0,145 0,01 0,123 2,50 a8 5,5 10,5 1 7,857 1,79 a9 0,55 1,05 0,1 0,843 5,36 a10 5,5 10,5 1 8,071 0,89 a11 0,145 0,195 0,01 0,172 1,18 На рис. 2.13 представлена динамика донорских антител. Величины скачков функции определяют объём вводимых в соответствующий момент времени донорских антител. Рис. 2.13. Динамика донорских антител 78 Результаты расчётов показывают, что использование рассмотренной методики позволяет строить управление в условиях неопределённости, когда значения пaраметров неизвестны, а их оценка корректируется по мере поступления новых экспериментальных значений. Набор параметров, который использовался при расчётах, не описывает какое-то конкретное заболевание, а характеризует типичную схему иммунной реакции на вторжение антигена. Под инфекционным заболеванием понимается отображение взаимоотношений, которые установились между антигеном и организмом. При этом организм с помощью иммунной системы способен оказывать противодействие патогенному влиянию антигена. Далее рассмотрим реализацию предложенных алгоритмов на примере конкретных заболеваний. 2.5. Моделирование управления иммунным ответом на основе реальных клинических данных Рассмотрим реализацию предложенного алгоритма в случае, когда значения параметров качественно воспроизводят динамику иммунной защиты oрганизма при острой пневмонии. После приведения параметров, взятых из работы [34], к безразмерному виду получим следующие значения: a1 = 0,25; a2 = 3,32⋅10−5 ; a3 = 7000; a4 = 0,05; a5 = 0,5; a6 = 10; a7 = 0,4; a8 = 1,7⋅10−3 ; τ = 0,5; m * = 0,1; v0 = 10−6 . В опорном решении a3 = 12000; τ = 0. В табл. 2.7 представлены результаты идентификации параметров модели после проведения 104 статистических испытаний. Средняя погрешность оценки параметров составляет 25,54%. Для сравнения в таблице приведены модельные значения, которые использовались ранее. Существенно отличаются значения параметров a1, a2, a8, тогда как величины параметров a3, a4, a5, a6, a7 примерно совпадают. При этом параметры a1, a2 характеризуют свойства антигена, поэтому при разных заболеваниях они могут существенно отличаться. 79 Таблица 2.7 Результаты идентификации параметров при пневмонии Параметр − i a + i a i h Модельные значения Оценка Погрешность, % a1 0,1 0,6 0,1 2 0,2 20,00 a2 10−5 6⋅10−5 10−5 0,8 5⋅10 −5 50,60 a3 6550 7550 100 10000 7350 5,00 a4 0,025 0,075 0,01 0,17 0,075 50,00 a5 0,25 0,75 0,1 0,5 0,55 10,00 a6 7,5 12,5 1 10 11,5 15,00 a7 0,15 0,65 0,1 0,12 0,35 12,50 a8 0,0005 0,003 0,0005 8 0,001 41,18 На рис. 2.14 представлено изменение вида управляющей функции в зависимости от момента времени начала лечения. а) б) Рис. 2.14. Управление иммунным ответом при пневмонии: а) 33, tN1−1 = б) 1 34 1 tN − = 80 Как видно из рисунка, вид управляющей функции качественно не изменился. Продолжительность лечения также лежит в пределах четырёх суток. Однако начинается иммунотерапия значительно позже. Это связано с более продолжительным доклиническим периодом заболевания. При этом существенно увеличился объём вводимых донорских антител. Далее рассмотрим идентификацию параметров на основе реальных клинических данных по вирусному гепатиту B [108]. У пациентов еженедельно измерялось количество вирусов. В табл. 2.8 представлены результаты идентификации параметров при острой форме заболевания. Проводилось 104 статистических испытаний. В данном случае значения параметров неизвестны, поэтому в таблице не представлена погрешность оценки, при этом наилучшая оценка определялась исходя из максимального согласования с клиническими данными. Таблица 2.8 Результаты идентификации параметров при гепатите B Параметр − i a + i a i h Оценка a1 0,1 1,1 0,1 0,9 a2 0,1⋅10−4 4,1⋅10−4 0,4⋅10−4 3,7⋅10−4 a3 550 10550 500 5050 a4 0,025 0,075 0,005 0,070 a5 0,25 0,75 0,1 0,75 a6 0,5 4,5 0,4 1,3 a7 0,05 0,55 0,05 0,40 a8 0,0005 0,3005 0,06 0,0305 На рис. 2.15 отмечены данные по динамике вирусов гепатита B, а также траектория, соответствующая полученной оценке параметров. Как видно из рисунка, фазовая траектория проходит вблизи экспериментальных значений. 81 Рис. 2.15. Динамика антигенов и клинические данные В табл. 2.9 представлено сравнение количественных характеристик при различных заболеваниях. Реальные заболевания от модельного существенно отличаются степенью повреждения органа-мишени и средней скоростью повреждения организма. При этом время выздоровления при гепатите значительно не отличается от модельного значения. Таким образом, реальные заболевания характеризуются значительным повреждением органа, что свидетельствует о необходимости лечения. Таблица 2.9 Количественные характеристики при различных заболеваниях Заболевание mmax I Tδ Модельное 0,025 0,06⋅10−3 21,59 Пневмония 0,260 1,50⋅10−3 50,56 Гепатит B 0,310 7,08⋅10−3 24,98 С помощью полученной оценки параметров построим управление. Опорное решение также будем определять исходя из минимизации функционала энергетической цены иммунного ответа. Вид управляющей функции 82 при вирусном гепатите B представлен на рис. 2.16. Опорное решение определялось при τ = 0 и a3 = 20050. Программа лечения заключается в непрерывном введении донорских антител в течение двух недель. Рис. 2.16. Управление иммунным ответом при гепатите Базовая математическая модель инфекционного заболевания содержит параметр m * , который характеризует максимальную степень повреждения организма, при которой ещё возможна нормальная работа иммунной системы. Принято m * = 0,1. Это значит, что если разрушено до 10% клеток органа, то работоспособность иммунологических органов не снижается. После превышения этого значения снижение работоспособности иммунной системы описано линейной функцией. Данный способ выглядит несколько искусственным. Поэтому зададим функцию ξ(m) следующим образом: ξ( ) 1 . 2 m = − m (2.49) В этом случае ухудшение работы иммунной системы описано квадратичной функцией вне зависимости от доли разрушенных антигенами клеток органа-мишени. В табл. 2.10 представлены результаты идентификации параметров при ξ(m) = 1 − m 2 . 83 Таблица 2.10 Результаты идентификации параметров при гепатите B (ξ(m) = 1 − m 2 ) Параметр − i a + i a i h Оценка a1 0,1 1,1 0,1 0,7 a2 0,1⋅10−4 4,1⋅10−4 0,4⋅10−4 2,1⋅10−4 a3 550 10550 500 3050 a4 0,025 0,075 0,005 0,030 a5 0,25 0,75 0,1 0,55 a6 0,5 4,5 0,4 0,9 a7 0,05 0,55 0,05 0,40 a8 0,0005 0,3005 0,06 0,1805 Таким образом, результаты, полученные для конкретного заболевания, хорошо согласуются с модельными расчётами, что свидетельствует о возможности применения алгоритма для прогноза течения заболевания и обоснования рекомендаций по выбору лечения. 2.6. Управление иммунным ответом при хронической и летальной форме заболевания Выше были представлены результаты для острой формы заболевания с выздоровлением. Управление в этом случае направлено на ускорение выздоровления. Рассмотрим применение описанной методики в случае хронической и летальной формы заболевания. Если при заражении здорового организма реализуется хроническая форма заболевания, то в этом случае в качестве опорного решения следует брать траекторию, соответствующую острой форме с выздоровлением. Для решения задачи (2.23) – (2.24) были использованы параметры, представленные в табл. 1.4 (хроническая форма). При хронической форме изменяются значения параметров, характеризующих скорость размножения антигенов и стимуляцию иммунной системы. Значения данных параметров при хронической форме заболевания меньше, чем при острой форме. Таким образом, в случае 84 хронической формы антигены обладают вялой динамикой, а стимуляция иммунной системы меньше, чем при острой форме. Реализация алгоритма при хронической форме заболевания представлена на рис. 2.17. Сплошными линиями показана управляющая функция. Момент времени начала лечения 1 33 1 tN − = сут. Опорное решение, которое соответствует острой форме заболевания с выздоровлением, изображено штриховой кривой. Множество значений (2.17) задавалось с интервалом времени Δt = 2 сут. Также на рисунке представлена динамика антигенов без управления (штрих-пунктир). Рис. 2.17. Управление при хронической форме заболевания Далее рассмотрим управление при хронической форме в условиях неопределённости. Результаты идентификации пaраметров в данном случае представлены в табл. 2.11. Наиболее точная оценка получена для параметров a2 и a3, наибольшая погрешность отмечается в оценке параметров a6 и a8. Средняя погрешность оценки параметров составляет 8,37%. Заметим, что при острой форме заболевания средняя погрешность оценки параметров составляет 8,06%. 85 Таблица 2.11 Результаты идентификации параметров (управление в условиях неопределённости при хронической форме) Параметр − i a + i a i h Оценка Погрешность, % a1 0,875 1,125 0,1 0,925 7,50 a2 0,55 1,05 0,1 0,800 0,00 a3 955 1050 10 985 1,50 a4 0,145 0,195 0,01 0,180 5,88 a5 0,25 0,75 0,1 0,550 10,00 a6 7,5 12,5 1 8,500 15,00 a7 0,095 0,145 0,01 0,130 8,33 a8 5,5 10,5 1 9,500 18,75 Управление в условиях неопределённости изображено на рис. 2.18 сплошными линиями. Для сравнения штриховыми линиями показана управляющая функция, построенная с заданными значениями параметров. Рис. 2.18. Управление при острой форме с возможным летальным исходом Острая форма заболевания с возможным летальным исходом характеризуется экспоненциальным ростом антигенов. Здесь управление должно оста- 86 новить рост антигенов и перевести форму заболевания в острую с выздоровлением. Интегрирование задачи (2.24) – (2.25) проводилось при значениях параметров, представленных в табл. 1.4 (летальный исход). На рис. 2.19 сплошными линиями изображена управляющая функция. Штриховой кривой показано опорное решение, соответствующее острой форме заболевания с выздоровлением. Также представлена динамика антигенов при реализации данного управления (штрих-пунктир). Рис. 2.19. Изменение концентрации антигенов v и донорских антител u со временем при острой форме с возможным летальным исходом Ввиду большого запаздывания при острой форме заболевания с возможным летальным исходом фазовая траектория концентрации антигенов на начальной стадии заболевания лежит ниже опорного решения, поэтому при рассмотренном способе управления Vi * > V(0, ti), i ∈ {1,2,…,6}. Вследствие этого введение донорских антител начинается при t = 6. Иммунотерапия позволяет переломить ход заболевания и приводит к выздоровлению. На рис. 2.20 представлено сравнение управления в условиях неопределённости (сплошные линии) с управляющей функцией, построенной c заданным набором значений пaраметров (штриховые линии). Отклонение управ- 87 ления в условиях неопределённости от управления при известных значениях параметров лежит в пределах 14%. Рис. 2.20. Управление при острой форме с возможным летальным исходом В табл. 2.12 представлены результаты идентификации пaраметров. Средняя погрешность оценки пaраметров составляет 5,73%. Таблица 2.12 Результаты идентификации параметров (управление в условиях неопределённости при острой форме с возможным летальным исходом) Параметр − i a + i a i h Оценка Погрешность, % a1 1,25 1,75 0,1 1,550 0,65 a2 0,55 1,05 0,1 0,650 15,58 a3 450 1450 100 879 0,11 a4 0,125 0,175 0,01 0,151 0,67 a5 0,25 0,75 0,1 0,393 21,40 a6 9,5 14,5 1 12,071 0,59 a7 0,095 0,145 0,01 0,115 4,17 a8 5,5 10,5 1 7,786 2,68 88 Таким образом, предложенный алгоритм позволяет строить управление при всех основных формах заболевания. 2.7. Дискретное управление в математической модели противовирусного иммунного ответа Базовая модель инфекционного заболевания описывает иммунный ответ достаточно схематично, поэтому рассмотрим применение описанной методики в случае более сложной модели. Будем использовать модель противовирусного иммунного ответа, которая является обобщением базовой модели. Математическая модель противовирусного иммунного ответа гуморального типа с учётом управления имеет вид . ξ( ) ( ) ( ) ( ) (1 ) , ξ( ) ( ) ( ) ( ) (1 ) , (1 ) , ξ( ) ( ) ( ) , , (1 ) , (1 ) , 3 2 3 3 3 4 2 9 3 0 3 0 3 0 3 1 2 5 2 6 2 6 2 6 2 7 2 8 1 7 1 8 1 4 1 5 1 5 1 6 6 7 3 7 3 8 3 5 3 7 1 3 4 5 a P a F V a F u dt dF a m M t a H t a B t a a P dt dP a m M t a H t a B t a a M H B a B dt dB a M H B a H a m H t a M t a a H M dt dH a V a M dt dM a C a m dt dm a V C m a C dt dC a C a V F a V a C m V dt dV f V B V B V B V B B B V B V B f V V V f V V V V f f V f f = − − + = − − − + − = − − − − + − − + − = − − − − = − = − = − − − = − − − − − (1.30) Начальные условия в момент инфицирования задаются следующим образом: ( ) ( ) ( ) 0, [ τ,0), τ max( , ). ( ) ( ) 0, [ ,0) , 0, 1, 2 6 3 0 1 5 0 M t H t B t t a a H t M t t a V v M C m H B P F V B B V f V V B = ∀ ∈ − = = ∀ ∈ − = = = = = = = = (2.49) 89 Опорное решение также будем формировать исходя из минимизации функционала энергетической цены иммунного ответа, который в модели противовирусного иммунного ответа имеет вид ( ) ( ) ( ). 0 = 2 9 ∫ − 2 6 − 2 6 − 2 6 T E a MV t a HB t a B t a (2.49) Результаты вычислительных экспериментов с управлением в математической модели противовирусного иммунного ответа опубликованы в работе [87]. На рис. 2.21 изображена зависимость энергетической цены иммунного ответа от коэффициента стимуляции иммунной системы при нулевом запаздывании. Опорное решение будем определять при a29 = 3,5 и τ = 0. Рис. 2.21. Энергетическая цена иммунного ответа в модели противовирусного иммунного ответа На рис. 2.22а) сплошными линиями показан вид управляющей функции. Также на рисунке представлена динамика антигенов при естественном течении заболевания (штриховая кривая) и в случае рассмотренной программы лечения (штрих-пунктир). Как видно из рисунка, через некоторое время после основного лечения необходимо повторное введение донорских антител для предотвращения рецидива заболевания. Если введение донорских антител заканчивается в момент времени t = 120, то наступает рецидив 90 заболевания (рис. 2.22б)). В этом случае концентрация антигенов при лечении (штрих-пунктир) не выходит на нулевое значение. а) б) Рис. 2.22. Управление иммунным ответом в модели противовирусного иммунного ответа Модель противовирусного иммунного ответа гуморального типа содержит двадцать четыре параметра. Для идентификации выбрано десять наиболее значимых параметров. Результаты идентификации параметров представлены в табл. 2.13. Максимальная погрешность оценки параметров составляет 10,0%, а средняя – 6,4%. 91 Таблица 2.13 Результаты идентификации параметров (математическая модель противовирусного иммунного ответа) Параметр − i a + i a i h Оценка Погрешность, % a1 0,085 1,115 0,006 0,109 9,0 a2 0,085 1,115 0,006 0,103 3,0 a5 0,045 0,055 0,002 0,055 10,0 a7 0,0095 0,0105 0,0002 0,0097 10,0 a12 0,75 0,85 0,02 0,750 6,7 a16 0,45 0,55 0,02 0,550 10,0 a19 0,165 0,175 0,002 0,171 0,6 a20 0,15 0,25 0,02 0,190 5,0 a22 0,35 0,45 0,02 0,370 7,5 a23 0,0045 0,0055 0,0002 0,0049 2,0 Таким образом, разработанная методика построения управления применима и к другим моделям инфекционных заболеваний. Качественный вид управляющей функции в базовой математической модели инфекционного заболевания совпадает с видом управления, построенного для модели противовирусного иммунного ответа. Выводы по главе 2 Сформулирован критерий управления иммунным ответом при острой форме инфекционного заболевания. Построена математическая модель для описания влияния иммунотерапии на динамику иммунного ответа. Предложен алгоритм дискретного управления иммунным ответом, с помощью которого построены и проанализированы программы лечения острой формы инфекционных заболеваний. 92 Сформулирован алгоритм идентификации пaраметров математических моделей иммунного ответа. С помощью предложенного алгоритма проведена идентификация пaраметров базовой модели инфекционного заболевания. Предложен алгоритм, позволяющий одновременно идентифицировать параметры и строить управление в математических моделях иммунного ответа при инфекционных заболеваниях. Предложенные алгоритмы протестированы на основе реальных клинических данных по динамике вирусного гепатита B. С помощью разработанной методики построено управление в математической модели противовирусного иммунного ответа. 93 3. МОДЕЛИРОВАНИЕ ЛЕЧЕНИЯ ХРОНИЧЕСКОЙ ФОРМЫ ЗАБОЛЕВАНИЯ 3.1. Общие свойства хронической формы заболевания Для хронической формы характерна вялая по сравнению с острой формой динамика антигенов. После заражения их количество медленно растёт, достигая максимума, далее также медленно снижается до минимума. Процесс роста и спада повторяется до установления в организме стационарного уровня антигенов V2. На рис. 3.1 показана динамика иммунного ответа при хронической форме заболевания. Как видно из рисунка, концентрация плазматических клеток и антител, а также доля разрушенных клеток органа тоже выходят на стационарный уровень. Решение задачи (1.26) – (1.27) проводилось при параметрах, представленных в табл. 1.4 (хроническая форма). Стационарное решение (1.18), характеризующее хроническую форму заболевания, после приведения модели к безразмерному виду запишется следующим образом: . ( ) ( ) , , ( ) , ( ) ( ) 1 7 3 4 5 8 4 5 6 1 2 2 2 1 2 2 3 4 5 8 1 3 4 2 5 8 2 1 3 4 5 8 4 5 1 2 2 a a a a a a a a a a a m a a f a a a a a a a a a a a s a a a a a a a a a v − − = = − − = − − = (3.1) В рамках линейной теории рассмотрим устойчивость к малым возмущениям, положив ε , ε , ε , ε . = * + 1 = * + 2 = * + 3 m = m* + 4 v v s s f f (3.2) Подставив выражения (3.2) в уравнения модели, получим ( ε ) ( ε ε ε ε ) ( )ε ε , ε ( ε ) ( ε ) ( ε ) 1 1 2 1 3 1 3 1 2 1 2 3 1 1 1 2 3 1 * * * * * * * * * * = + − + + + ≈ − − = + − + + = a v a f v f v a a f a v & a v a f v ε ( τ) ε ε ( τ), ( ε ( τ) ε ( τ) ε ε ) ( ε 1) ε ( ε ( ) ) ( ε ( τ) ) ( ε 1) 3 1 5 2 3 3 3 1 3 1 3 5 2 2 3 3 1 5 2 ≈ − − + − = + − + − + − + − ≈ = + − + − − + − = * * * * * * * * * * a f t a a v t a f v f t v t a s & a f t τ v t a s ( ε ε ε ε ) ε ε ( )ε , ε ( ε ε ) ( ε ) ( ε ) ( ε ε ) 8 1 3 1 3 8 1 4 2 4 8 3 3 4 2 3 8 3 1 4 2 3 * * * * * * * * * * * * − + + + ≈ − + − + = + − − − + + = + − − − a f v f v a f a a a v & a s f a f v a s f ε ( ε ) ( ε ) ε ε . 4 6 1 a7 m 4 a6 1 a7 4 & = a v* + − * + ≈ − (3.3) 94 Таким образом, малые возмущения удовлетворяют системе ε, ε ρ ρ A dt d = (3.4) где . 0 0 ( ) 0 0 0 0 6 7 8 4 4 8 λ τ 5 3 λ τ 3 1 2 2 ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − − − + − − − = * * − * − * * * a a a f a a a v a f e a a v e a a f a v A (3.5) Необходимо проанализировать значения корней характеристического уравнения A− λE = 0. (3.6) Характеристический квазиполином имеет вид ( λ)( λ λ λ ( λ ) ) 0, 3 2 λ τ − 7 + − − − + + − = − a a b d g f e (3.7) где , 4 5 8 2 a = a + a + a v ( ) , 5 8 2 4 1 8 2 b = a a v + a − a a v , 1 5 8 2 d = a a a v , 3 4 2 g = a a v . 1 3 4 2 f = a a a v (3.8) Достаточным условием асимптотической устойчивости стационарного решения (3.1) при a5τ ≤ 1 является выполнение условий [6]: , a3a4 > a5a8 , a1 > a2 . τ ( ) ( ) ( ) τ( ) 0 1 3 4 5 8 8 3 4 1 4 1 2 2 4 5 1 1 5 3 4 5 8 5 1 2 2 4 1 3 4 2 5 8 3 4 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − + − − < ⎥ < ⎦ ⎤ ⎢ ⎣ ⎡ − − − − < + − a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a (3.9) 95 а) б) в) г) Рис. 3.1. Динамика иммунного ответа при хронической форме заболевания: а) антиген, б) плазматические клетки, в) антитела, г) доля разрушенных клеток 96 На рис. 3.2 показана область устойчивости стационарного решения (3.1) в зависимости от параметров a3, a4. Данные параметры отвечают за скорость производства плазматических клеток и антител и наиболее сильно влияют на характер течения заболевания. Если состояние хронической формы неустойчиво, то в модели реализуется острая форма заболевания с выздоровлением или с возможным летальным исходом. Область устойчивости расположена над изображённой кривой, полученной из первого соотношения (3.9). Рис. 3.2. Область устойчивости хронической формы Моделирование инфекционного заболевания как конфликта между иммунной системой и размножающимся антигеном показало, что хронические формы связаны с вялым иммунным ответом. Оказалось, что значительное увеличение дозы инфицирования может привести к переводу хронической формы в острую форму с выздоровлением. При анализе базовой математической модели инфекционного заболевания был сделан вывод о том, что хронические инфекции можно лечить обострением. Это привело к построению модели протекания иммунного ответа на проникновение в oрганизм двух различных антигенов [8], или биинфекции. На основе этой модели в работе [7] предложено управление иммунной реакцией при хронической форме за- 97 болевания. Лечение хронической формы здесь основано на взаимодействии двух антигенов: возбудителя болезни и специально вводимого биостимулятора. Это даёт возможность эффективно стимулировать иммунный ответ. После прекращения инъекций биостимуляторов образуется большое количество антител, в результате концентрация антигенов резко падает, вследствие чего наступает выздоровление. В работе [17] на основе данного подхода с использованием принципа максимума Понтрягина определяются оптимальные значения биостимуляторов. Сформулируем ряд других способов управления иммунной реакцией при хронической форме заболевания. 3.2. Управление иммунным ответом на основе анализа стационарного решения Вначале сформулируем цель управления иммунным ответом при хронической форме заболевания: ( ) 0, (0) . 2 v t v v t ⎯ ⎯→ = →∞ (3.10) Таким образом, систему необходимо перевести из состояния устойчивой хронической формы в состояние здорового организма. Для решения поставленной задачи рассмотрим стационарное решение базовой модели инфекционного заболевания с управлением: 0. ( ) 0, ξ( ) ( τ) ( τ)θ( τ) ( 1) 0, 0, 6 7 4 8 3 5 1 2 − = − − + = − − − − − = − = a v a m a s f a f v u a m f t v t t a s a v a f v (3.11) Из первого уравнения получим ( ) 0, v a1 − a2 f = (3.12) откуда следует соотношение для концентрации антител . 2 1 a a f = (3.13) 98 Учитывая, что управление должно обеспечить равенство концентрации антигенов нулю (v = 0), из второго уравнения получим s = 1. Из третьего уравнения выведем соотношение для управляющей функции: ( ) 1 . 2 1 4 4 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − = − * a a u a f s a (3.14) При значениях параметров, представленных в табл. 1.4 (хроническая форма), получим u * = 0,0425. На рис. 3.3 представлен вид управляющей функции, соответствующий данному способу управления. Программа лечения заключается в непрерывном введении донорских антител в течение 25 суток. Введение донорских антител прекращается, когда их концентрация опускается ниже некоторой малой величины ε, в данном случае ε = 10−9 . Рис. 3.3. Вид управляющей функции Таким образом, зная значения параметров модели, можно найти величину управления. Однако в реальных условиях параметры системы неизвестны, поэтому необходимо использовать методы, позволяющие строить управление в условиях неопределённости. 99 3.3. Дискретное управление иммунным ответом при хронической форме заболевания Рассмотрим построение управления при хронической форме на основе подхода, изложенного в главе 2. Запишем управляемую модель иммунного ответа при инфекционном заболевании в общем виде: ( ) , [ τ,0], ( , , ( τ), ), [0, ], = хр. ∈ − = − ∈ t t t u t T dt d x x F α x x x (3.15) где xхр. соответствует состоянию устойчивой хронической формы. Будем считать, что экспериментальные значения можно получить в узлах сетки : , 0, , . ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ Π = = Δ = Δ = N T t t i t i N t i i (3.16) Алгоритм, опубликованный в работах [68, 86], заключается в следующем. Сначала задаются величины допустимого отклонения расчётных значений от экспериментальных данных εj , j = 1, …, m. Затем случайным образом задаётся K наборов значений параметров. При t ∈ Π определяются допустимые наборы пaраметров, то есть удовлетворяющие критерию ( , ) ( ) ε , 0, , 1, , 1, , ( ) эксп. y t yj t i j i N k Ki j m k j i α − < = = = (3.17) где Ki – количество наборов пaраметров в момент ti . Таким образом, после измерения значений фазовых переменных при t = ti необходимо найти Ki решений задачи (3.15) на отрезке [0, ti], i = 0, …, N, и определить допустимые наборы пaраметров. В качестве оценки параметров при t = ti выбирается среднее значение допустимых наборов: , 1, , 0, , 1 ( ) ( ) j L i N J a a i J k k j i j i = = = ∑ = (3.18) где Ji – количество допустимых наборов значений пaраметров в момент времени ti ; Ji ≤ Ki , i = 0, …, N; Ki = Ji − 1, i = 1, …, N; K0 = K , где K – первоначаль- 100 ное количество наборов пaраметров; Ji = Ki − Hi , i = 0, …, N, где Hi – количество неприемлемых наборов пaраметров в точке ti . В качестве опорного решения будем задавать убывающую часть траектории из решения задачи со значениями пaраметров, соответствующими острой форме заболевания: , 1, , , 1, 1, . 1 N V1 V2 V i N V V i i = i ≥ i = − ≤ * * + * * (3.19) Для каждого допустимого набора пaраметров подбирается значение управления i k ui , i 1,N, k 1,K ( ) = = . Управляющая функция представляет собой среднее значение найденных величин: , 1, . 1 ( ) i N J u u i J k k i i i = = ∑ = (3.20) Таким образом, предложенный алгоритм позволяет строить управление иммунной реакцией при хронической форме заболевания в условиях неопределённости. Рис. 3.4. Динамика антигенов Рассмотрим реализацию алгоритма. На рис. 3.4 представлена динамика антигенов при опорном решении (сплошная кривая, также на графике отме- 101 чены опорные точки (3.19) при Δt = 0,5 сут., которые используются для построения управления), в начальной стадии при острой форме заболевания (штриховая кривая), а также в случае устойчивой хронической формы (штрих-пунктир). На рис. 3.5 показан вид управляющей функции. Сплошными линиями изображена управляющая функция, построенная в условиях неопределённости (при одновременной идентификации параметров), а штриховыми – с заданным набором значений пaраметров. В данном случае максимальное отклонение управления в условиях неопределённости от управляющей функции, построенной с заданным набором значений параметров, лежит в пределах 41%. Однако отклонение объёма введения донорских антител не превышает 7%. Рис. 3.5. Управляющая функция Результаты идентификации пaраметров при K = 104 представлены в табл. 3.1. Средняя погрешность оценки пaраметров составляет 6,17%, максимальная погрешность – 15%. 102 Таблица 3.1 Результаты идентификации параметров при хронической форме Параметр − i a + i a i h Оценка Погрешность, % a1 0,875 1,125 0,1 1,037 3,70 a2 0,55 1,05 0,1 0,800 0,00 a3 955 1055 10 1023 2,30 a4 0,145 0,195 0,01 0,175 2,94 a5 0,25 0,75 0,1 0,425 15,00 a6 7,5 12,5 1 11,500 15,00 a7 0,095 0,145 0,01 0,125 4,17 a8 5,5 10,5 1 7,500 6,25 На рис. 3.6 представлен вид управляющей функции в зависимости от продолжительности промежутков введения донорских антител в случае непрерывного введения. Как видно из рисунка, при уменьшении интервалов времени, в течение которых вводятся донорские антитела, вид управляющей функции не изменяется, однако снижается продолжительность лечения (наблюдается эффект праймирования). На рис. 3.7 представлено изменение концентрации донорских антител при дискретном введении в зависимости от продолжительности интервалов времени между соседними инъекциями. Как видно из рисунка, при уменьшении интервала времени между инъекциями наблюдается увеличение числа инъекций. Сравнение полученных программ лечения будем проводить по средней скорости повреждения организма, объёму введения донорских антител, времени, необходимому для полного выздоровления, и продолжительности лечения. 103 а) б) г) Рис. 3.6. Вид управляющей функции: а) Δt = 1 сут., б) Δt = 0,5 сут., в) Δt = 0,25 сут. 104 В табл. 3.2 представлены значения выбранных характеристик для рассмотренных программ лечения. Таблица 3.2 Сравнительные данные программ лечения Программа лечения I Q Tδ Tл Рис. 3.4 Δt = 1 сут. 9,74⋅10−5 9,47 5,80 2 Δt = 0,5 сут. 5,66⋅10−5 7,63 5,24 1 Δt = 0,25 сут. 3,58⋅10−5 6,08 4,96 0,5 Рис. 3.6 Δt = 1 сут. 8,45⋅10−5 5,06 5,62 1 Δt = 0,5 сут. 5,39⋅10−5 4,23 5,21 0,5 Δt = 0,25 сут. 3,67⋅10−5 4,65 4,97 0,5 Из анализа представленных результатов можно сделать вывод о том, что наиболее приемлемым является вариант с продолжительностью введения донорских антител, равной четверти суток, так как в этом случае наблюдается наименьший объём введения донорских антител, наименьшая скорость повреждения организма, наименьшая продолжительность лечения, а также раньше наступает выздоровление. При инъекциях лекарственных препаратов наименьший объём введения донорских антител наблюдается при лечении с интервалом времени между инъекциями, равным половине суток. Программа лечения с интервалом времени между соседними инъекциями, равным четверти суток, характеризуется быстрейшим выздоровлением и наименьшей скоростью повреждения организма. Наименьшая продолжительность лечения имеет место при Δt = 0,5 сут. и Δt = 0,25 сут. 105 а) б) в) Рис. 3.7. Динамика донорских антител: а) Δt = 1 сут., б) Δt = 0,5 сут., в) Δt = 0,25 сут. 106 Отсюда можно сделать вывод o том, что более эффективным является способ управления с интервалом времени между инъекциями, равным четверти суток. 3.4. Моделирование лечения хронической формы заболевания с помощью алгоритмов управления с обратной связью Рассмотренный способ идентификации параметров и управления при хронической форме заболевания даёт достаточно большую погрешность при построении управления. Это может быть связано с малой информативностью фазовых траекторий при хронической форме заболевания. В отличие от острой формы, для которой характерна выраженная динамика иммунного ответа, хроническая форма обладает вялой динамикой. Поэтому в этом случае воспользуемся алгоритмом управления с обратной связью. Сначала рассмотрим управление с обратной связью на примере системы «хищник-жертва», поскольку она является типичным представителем класса моделей популяционной динамики, имеет достаточно простую структуру и небольшое число пaраметров, что позволяет провести её детальное аналитическое исследование. Кроме того, система «хищник-жертва» применялась для моделирования иммунного ответа [93]. 3.4.1. Управление с обратной связью в системе «хищник-жертва» Задача управления для классической системы типа «хищник-жертва» может быть записана следующим образом [32]: α β ( , , ), X & = 1X − 1XY +U t X Y α β , 0, Y & = − 2Y + k 1XY t > (3.21) 0 1 α ,β ,α 0, < k < , 1 1 2 > (0) , (0) , , 0. X = X0 Y =Y0 X0 Y0 > (3.22) Будем полагать, чтo значения пaраметров задачи (3.21) – (3.22) неизвестны, но значения функций X(t) и Y(t) могут быть получены (измерены) в 107 любой момент времени c необходимой точностью (кибернетическая модель «белый ящик»), а управление U(t,X,Y) означает возможность добавлять (убавлять) в систему извне любое необходимое количество «жертвы». Цель управления заключается в обеспечении асимптотической стабилизации функции Y(t) около заданного значения: ( )⎯ ⎯→ = > 0. →∞ Y t A const t (3.23) Приведём задачу (3.21) – (3.22) к безразмерному виду, выбрав в качестве характерного значения «хищника» и «жертвы» величину A, а характерного времени – величину 1 α2 − Δt = , сохраняя прежнее обозначение для времени, получим x & = αx −βxy + u(t, x, y), y & = −y + kβxy, (3.24) 0 < k <1, α,β > 0, (0) , (0) , , 0. x = x0 y = y0 x0 y0 > (3.25) ( )⎯ ⎯→1, t→∞ y t (3.26) где . α ( , , ) , ( , , ) α β , β α α α 2 2 1 2 1 A U t X Y = = A u t x y = (3.27) Сначала проанализируем вариант u(t,x,y) = const с помощью первого метода Ляпунова. В этом случае нетривиальным стационарным решением задачи (3.24) – (3.25) является: , β α 1 1 , 1, β 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = = = − * * * k y u k x (3.28) и корни соответствующего характеристического уравнения имеют вид β. 2 α β , 2 α β λ 2 1,2 ⎟ − ⎠ ⎞ ⎜ ⎝ ⎛ − ± = − = D D (3.29) При α < β это решение является устойчивым узлом в случае D > 0 и устойчивым фокусом при D < 0. В случае α > β имеем соответственно неустой- 108 чивый узел и неустойчивый фокус. Таким образом, в этом варианте устойчивое решение получится только в случае выполнения условия α < β. Поскольку параметры системы (3.21) неизвестны, достижение цели управления (3.23) возможно только при использовании алгоритмов управления с обратной связью. Для реализации такого алгоритма будем использовать ПИД-регулятор (пропорциональный, интегрпальный, дифференциальный), который широко используется при управлении различными техническими системами [4]: ( , , ) ( ) (τ) τ ( ), , , , Д П И Д 0 П И u t x y K e y K e d K e y K K K const t = ⋅ + ⋅ ∫ + ⋅ & = (3.30) где e(y) = 1 − y(t) – функция рассогласования. Полученные результаты опубликованы в статье [67]. Для исследования устойчивости системы (3.24) c управлением вида (3.30) методом Ляпунова предварительно преобразуем её, продифференцировав первое уравнение: x & = z, z & = (α −βy)z +βxy(1− kβx) + u &(t, x, y,z), y & = −y + kβxy, (3.31) где ( , , , ) (1 β )( (1 β )) (1 ) β . П Д И Д u & t x y z = −y − k x −K + K − k x + K − y − K k zy (3.32) Соответствующее характеристическое уравнение имеет вид λ λ λ 0, 2 3 2 1 3 + a + a + a = (3.33) где β(1 ) α, β(1 ), β. 1 Д 2 П 3 И a = + K k − a = + K k a = K k (3.34) Воспользуемся критерием Рауса-Гурвица, согласно которому все корни характеристического уравнения (3.33) имеют отрицательную действительную часть тогда и только тогда, когда выполняютс я следующие условия: β(1 ) α 0, T1 = a1 = + KД k − > (3.35) 109 (β(1 ) α)β(1 ) β 0. 1 1 2 3 Д П И 3 2 1 2 = = a a − a = + K k − + K k − K k > a a a T (3.36) Следствием соотношений (3.35) – (3.36) являются следующие ограничения: 1 , β 1 α Д Д * =⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ > − K k K (3.37) (β(1 ) α)(1 ) . 1 И Д П * < + − + = Ki K k K k k K (3.38) Таким образом, основное функциональное назначение составных частей (слагаемых) ПИД-регулятора заключается в следующем: 1. дифференциальная составляющая обеспечивает устойчивость управления в случае α > β (условие (3.37)); 2. интегральное слагаемое даёт необходимое значение величины управления * К ⋅ ∫ e d ⎯→⎯∞ →u t t 0 И (τ) τ (3.39) и должно удовлетворять условию (3.38); 3. пропорциональная составляющая позволяет регулировать скорость выхода нa стационарное решение. Таким образом, для рассматриваемой системы недостаточно ПИрегулятора и только наличие Д-компонеты обеспечивает возможность устойчивого управления для любого соотношения пaраметров системы. Полученные результаты можно использовать при идентификации параметров реальной системы типа «хищник-жертва», для которой имеется возможность проводить серию экспериментов. В случае отсутствия управления тип особой точки «центр» и размерное время, за которое система совершает полный цикл, определяется из соотношения . α α 2π 1 2 Tc = (3.40) 110 Поскольку эта величина достаточно просто определяется из экспериментальных данных, то можно получить следующее соотношение для параметров системы: . 2π α α 2 1 2 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = Tc (3.41) Кроме того, эксперимент даёт возможность оценить снизу и сверху отношение α1/β1. При применении стабилизирующего управления, обеспечивающего значение Y* = A, количество «жертвы» и интегральная составляющая ПИД-регулятора стремятся к значениям (3.28): , β α β 1 2 k k A X* = = (3.42) . β 1 α β α 1 1 1 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ * = − A k k A U (3.43) Таким образом, для определения четырёх пaраметров системы (3.21) α1, α2, β1, k имеем три соотношения (3.41) – (3.43). Если провести две серии экспериментов со стабилизацией численности «хищника» около значений Ai (i = 1, 2), получим , 1,2. β 1 α β α 1 1 ( ) 1 = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ * = − A i k k A U i i i (3.44) Тогда из уравнений (3.44) имеем , (1) (2) 1 2 * − * − = U U A A k (3.45) , β α (1) 1 1 1 * = A − k ⋅U (3.46) и система уравнений (3.41), (3.42), (3.46) легко разрешается относительно неизвестных α1, α2, β1. Таким образом, использование управления в виде ПИД-регулятора позволяет эффективно решать обратную задачу определения пaраметров классической системы «хищник-жертва». 111 3.4.2. Управление иммунным ответом при хронической форме заболевание с помощью ПИД-регулятора Выберем следующую цель управления: ( )⎯ ⎯→0. t→∞ v t (3.47) Таким образом, с помощью управления необходимо перевести систему в состояние здорового организма. Будем формировать управляющую функцию в виде ПИД-регулятора. В основном ПИД-регуляторы используются в тех случаях, когда необходимо перевести систему из одного стационарного состояния в другое. В случае управления иммунным ответом при хронической форме заболевания необходимо перевести систему из состояния устойчивой хронической формы в состояние здорового организма. В качестве функции рассогласования выберем динамику концентрации антигенов, тогда управление будет выглядеть следующим образом: ( ) ( ) , , , . Д П И Д 0 П И u t K v K v x dx K v K K K const t = + ∫ + & = (3.48) Для вычисления интегральной составляющей использовалась формула трапеций: ( ) , 2 1 ( ) 1 1 0 0 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∫ ≈ + + ∑ − = n i n i t v s ds h v v v (3.49) где vn = v(t). В рамках линейной теории рассмотрим устойчивость управляемой системы к малым возмущениям. Поскольку управляющая функция содержит интегральную составляющую, уравнение для концентрации антител необходимо продифференцировать второй раз. Пусть f = z & , тогда ξ( ) ( τ) ( τ) ( 1) , , , ξ( ) ( τ) ( τ) ( 1), , 2 3 4 4 5 8 1 8 2 8 6 7 3 5 1 2 z a a m f t v t a a s a zv a a v f a a f v u m a v a m f z s a m f t v t a s v a v a f v & & & & & & = − − − − − − + + = − = = − − − − = − (3.50) 112 где ( ) ( 2 ). 2 2 1 2 2 2 2 П 1 2 И Д 1 u & = K v a − a f + K v + K v a − a a f − a z − a f (3.51) Сделав замену переменных 1 2 3 4 5 v = ε , s =1+ ε , f =1+ ε , m = ε , z = ε (3.52) получим ε ε ( ( ) ( ( ) ) ε ε . ε ε ε , ε ε , ε ε ε , ε ( )ε , 5 1 3 4 1 2 П Д 1 2 И 8 2 4 5 5 4 4 6 1 7 4 3 5 2 3 1 5 2 1 1 2 1 a a a a K K a a K a a a a a a a a a a = + − + − + − − − = − = = − = − & & & & & (3.53) Характеристическое уравнение имеет вид λ( λ)( λ)( λ)( λ) 0. a1 − a2 − a4 + a5 + a7 + = (3.54) Корни характеристического уравнения λ 0, λ , λ , λ , λ . 1 = 2 = a1 − a2 3 = −a4 4 = −a5 5 = −a7 (3.55) Следовательно, для устойчивости необходимо выполнение условия a1 < a2. Параметры ПИД-регулятора не входят в корни характеристического уравнения, поэтому они не влияют на устойчивость системы. Следовательно, на значения параметров ПИД-регулятора ограничения не накладываются. Поэтому рассмотрим изменение вида управляющей функции при варьировании параметров ПИД-регулятора. На рис. 3.8 представлено влияние на управляющую функцию пропорциональной составляющей. В качестве параметров интегрального и дифференциального слагаемых были взяты следующие значения: KИ = 15, KД = 2. Сплошной кривой показана управляющая функция при KП = 2, а штриховой – при KП = 50. Как видно из рисунка, изменение пропорционального слагаемого слабо влияет на управление. 113 Рис. 3.8. Зависимость управления от пропорционального слагаемого Рис. 3.9 демонстрирует изменение вида управляющей функции при варьировании интегрального параметра (KП = KД = 2): KИ = 2 (сплошная кривая), KИ = 5 (пунктирная кривая), KИ = 12 (штриховая кривая), KИ = 20 (штрих-пунктир). Интегральная составляющая оказывает наиболее существенное влияние на вид управления, так как она даёт необходимое значение величины управления. Рис. 3.9. Зависимость управления от интегрального слагаемого 114 В результате вычислительных экспериментов установлено, что изменение параметра дифференциальной составляющей не влияет на вид управляющей функции. Таким образом, дифференциальную составляющую можно опустить. Поэтому далее будем использовать интегральный регулятор. Для выбора коэффициента интегрального регулятора рассмотрим следующую оптимизационную задачу: ( ) ( ) ( ) min. И 0 τ τ 3 К T T E = a ∫ v t f t dt + ∫ u t dt → − − (3.56) В результате вычислительных экспериментов установлено, что оптимальным является значение параметра KИ = 22. Зависимость величины E от параметра КИ изображена на рис. 3.10. Рис. 3.10. Зависимость E от параметра КИ На рис. 3.11 изображена динамика антигенов. Сплошной кривой показан полученный с помощью управления выход системы на состояние здорового организма, а штриховой – наличие устойчивой хронической формы. 115 Рис. 3.11. Динамика антигенов На рис. 3.12 показан вид управляющей функции при указанных значениях параметров. Рис. 3.12. Управляющая функция Таким образом, в рассмотренном варианте управления наблюдается выход на стационарное решение, соответствующее состоянию здорового организма. 116 3.5. Идентификация параметров модели при хронической форме заболевания Полученные результаты можно использовать для идентификации параметров модели. Учитывая, что значения фазовых переменных, соответствующие стационарному решению (3.1), можно измерить в ходе эксперимента, получаем четыре уравнения для определения параметров системы. Введём обозначение , , , . 4 7 6 3 8 4 2 3 5 1 2 1 b a a b a a b a a b a a = = = = (3.57) Здесь отношения параметров подобраны таким образом, что в числителе b1 и b4 стоит параметр, отвечающий за скорость роста соответствующей фазовой переменной, а в знаменателе – за скорость убывания. Подставив выражения (3.57) в стационарное решение (3.1), получим , , , . ( ) ( 1) 2 1 2 4 2 3 2 1 3 2 2 1 3 2 3 2 1 2 f b m b v b b b b b s b b b b b b v = = − − = − − = (3.58) Величину (3.14) также можно определить в ходе эксперимента. Таким образом, получаем систему уравнений ( 1). , , , , ( ) ( 1) 4 1 2 4 2 2 1 3 2 1 3 2 2 1 3 2 3 2 1 2 = − = = − − = − − = * u a b m b v f b b b b b b s b b b b b b v (3.59) В данной системе неизвестными являются параметры b1, b2, b3, b4, a4. Разрешив систему (3.59) относительно неизвестных b1, b2, b3, b4, a4, получим значения a1 /a2, a5 /a3, a4, a6 /a7, a8. Учитывая, что параметры, отвечающие за скорость убывания фазовых переменных, незначительно влияют на изменение решения, их значения можно задавать из теоретических соображений. Задав значения параметров a2, a5, a7 можно идентифицировать параметры a1, a3, a4, a6, a8. 117 Таким образом, использование управления в виде интегрального регулятора позволяет решать задачу идентификации параметров базовой модели инфекционного заболевания при хронической форме. Выводы по главе 3 Для хронической формы заболевания найдено управление, позволяющее при заданных значениях параметров перевести систему из состояния хронической формы в состояние здорового организма. Сформулирована задача дискретного управления иммунным ответом при хронической форме заболевания в условиях неопределённости, а также предложен алгоритм её решения, позволяющий строить управление одновременно с идентификацией параметров. Рассмотрено управление с обратной связью в виде интегрального регулятора, найдено оптимальное значение коэффициента интегрального регулятора, построена управляющая функция. 118 ЗАКЛЮЧЕНИЕ 1. На основе базовой модели инфекционного заболевания путём расширения пространства фазовых переменных и модификации уравнений построена математическая модель для описания влияния иммунотерапии на динамику иммунного ответа. 2. Сформулирован критерий управления заболеванием, который обеспечивает «идеальный» иммунный ответ, соответствующий минимальному расходу энергии на взаимодействие с антигеном. 3. Предложены алгоритмы, позволяющие строить управление с заданным набором значений параметров модели, а также в условиях неопределённости, когда значения параметров неизвестны, а их оценка корректируется по мере поступления новых клинических данных. Создан программный инструментарий для решения задачи управления иммунным ответом при инфекционных заболеваниях. С помощью программной реализации алгоритмов спрогнозирован ход течения заболевания, а также получены результаты по управлению иммунным ответом. 4. Построены программы лечения острой и хронической формы инфекционных заболеваний. Проведён анализ их эффективности. Предложенные алгоритмы протестированы на основе реальных клинических данных по динамике вирусного гепатита B. 119 СПИСОК ЛИТЕРАТУРЫ 1. Асаченков А.Л. Простейшая модель влияния температурной реакции на динамику иммунного ответа // Математическое моделирование в иммунологии и медицине / под ред. Г.И. Марчука. – Новосибирск: Наука, 1982. – С. 40 – 44. 2. Асаченков А.Л., Белых Л.Н. Исследование математической модели вирусного заболевания // Математические методы в клинической практике. – Новосибирск: Наука, 1978. – С. 19 – 26. 3. Асаченков А.Л., Марчук Г.И. Уточнённая математическая модель инфекционного заболевания // Математическое моделирование в иммунологии и медицине / под ред. Г.И. Марчука. – Новосибирск: Наука, 1982. – С. 44 – 59. 4. Афанасьев В.Н., Колмановский В.Б. Носов В.Р. Математическая теория конструирования систем управления. 3-е изд., испр. и доп. – М.: Высш. шк., 2003. – 614 с. 5. Бард Й. Нелинейное оценивание параметров. – М.: Статистика, 1979. – 352 с. 6. Белых Л.Н. Анализ математических моделей в иммунологии. – М.: Наука, 1988. – 192 с. 7. Белых Л.Н. Математическая модель биинфекции и лечение хронических форм болезни обострением в её рамках // Математическое моделирование в иммунологии и медицине / под ред. Г.И. Марчука. – Новосибирск: Наука, 1982. – С. 33 – 40. 8. Белых Л.Н. Математическая модель присоединённого заболевания // Математические модели заболеваний и методы обработки медицинской информации / под ред. Г.И. Марчука. – Новосибирск: Наука, 1979. – С. 32 – 37. 120 9. Белых Л.Н. О численном решении моделей заболеваний // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 291 – 297. 10. Белых Л.Н., Марчук Г.И. Качественный анализ простейшей математической модели инфекционного заболевания // Математическое моделирование в иммунологии и медицине / под ред. Г.И. Марчука. – Новосибирск: Наука, 1982. – С. 5 – 27. 11. Белых Л.Н., Марчук Г.И., Петров Р.В. О некоторых подходах к математическому моделированию в иммунологии // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 5 – 22. 12. Бернет Ф. Клеточная иммунология. – М.: Мир, 1971. – 542 с. 13. Болодурина И.П., Луговскова Ю.П. Математическая модель управления иммунной системой // Обозрение прикладной и промышленной математики. – 2008. – Т. 15. Вып. 6. – С. 1043 – 1044. 14. Болодурина И.П., Луговскова Ю.П. Моделирование оптимальных стратегий лечения с учётом энергетических затрат противоинфекционной защиты // Математическая теория систем: материалы международной конференции. – М.: ИСА РАН, 2009. – С. 142 – 149. 15. Болодурина И.П., Луговскова Ю.П. Моделирование стохастических колебаний нелинейной системы дифференциальных уравнений, описывающих динамику иммунного ответа при инфекционных заболеваниях // Вестник ОГУ. – 2010. – № 9(115). – С. 23 – 29. 16. Болодурина И.П., Луговскова Ю.П. Оптимальное управление динамикой взаимодействия иммунной системы человека с инфекционными заболеваниями // Вестник СамГУ. Естественнонаучная серия. – 2009. – № 8(74). – С. 138 – 153. 17. Болодурина И.П., Луговскова Ю.П. Оптимальное управление иммунным ответом при хронической форме инфекционных заболеваний // Обозре- 121 ние прикладной и промышленной математики. – 2009. – Т. 16. Вып. 2. – С. 296 – 297. 18. Болодурина И.П., Луговскова Ю.П. Оптимальное управление иммунологическими реакциями организма человека // Проблемы управления. – 2009. – № 5. – С. 44 – 52. 19. Бочаров Г.А. Исследование математической модели противовирусного T-клеточного иммунного ответа // Вычислительные процессы и системы. Вып. 3 / под ред. Г.И. Марчука. – М.: Наука, 1985. – С. 199 – 207. 20. Бочаров Г.А. Численные эксперименты с математической моделью Tклеточного иммунного ответа // Вычислительные процессы и системы. Вып. 3 / под ред. Г.И. Марчука. – М.: Наука, 1985. – С. 208 – 218. 21. Бочаров Г.А., Марчук Г.И. Прикладные проблемы математического моделирования в иммунологии // Ж. вычисл. матем. и матем. физ. – 2000. – Т. 40. № 12. – С. 1905 – 1920. 22. Гноенский Л.С., Каменский Г.А., Эльсгольц Л.Э. Математические основы теории управляемых систем. – М.: Наука, 1969. – 512 с. 23. Горбунов А.Д., Попов В.Н. О методах типа Адамса приближённого решения задачи Коши для обыкновенных дифференциальных уравнений с запаздыванием // Численные методы решения дифференциальных и интегральных уравнений и квадратурные формулы. – М.: Наука, 1964. – С. 135 – 148. 24. Дибров Б.Ф., Лифшиц М.А., Волькенштейн М.В. Математическая модель иммунной реакции I // Биофизика. – 1976. – Т. 21. – С. 905 – 909. 25. Дибров Б.Ф., Лифшиц М.А., Волькенштейн М.В. Математическая модель иммунной реакции II // Биофизика. – 1977. – Т. 22. – С. 313 – 317. 26. Дружченко В.Е. Математическая модель противовирусного иммунного ответа гуморального типа // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 168 – 174. 122 27. Зуев С.М. Математические модели заболеваний и анализ экспериментальных данных. – М.: Наука, 1987. – 192 с. 28. Зуев С.М. Определение параметров моделей иммунного ответа по данным наблюдений // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 298 – 308. 29. Зуев С.М. Статистическое оценивание параметров математических моделей заболеваний. – М.: Наука, 1988. – 176 с. 30. Зуев С.М., Дружченко В.Е. Идентификация простейшей математической модели вирусного заболевания // Вычислительные процессы и системы. Вып. 3 / под ред. Г.И. Марчука. – М.: Наука, 1985. – С. 219 – 227. 31. Иванов В.В., Яненко В.М., Фонталин Л.Н., Нестеренко В.Г. Моделирование идиотип-антиидиотипических взаимодействий иммунной сети с учётом деления лимфоцитов на субпопуляции // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М: Мир, 1986. – С. 123 – 135. 32. Ильичёв В.Г. Устойчивость, адаптация и управление в экологических системах. – М.: Физматлит, 2009. – 192 с. 33. Каркач А.С., Романюха А.А. Энергетический критерий качества иммунной защиты и патогенность микроорганизмов // Автоматика и телемеханика. – 2003. – № 6. – С. 162 – 172. 34. Каркач А.С., Санникова Т.Е., Романюха А.А. Математическое моделирование адаптации иммунной системы. Энергетическая цена приспособленности // Вычислительная математика и математическое моделирование: труды международной конференции / под ред. В.П. Дымникова. Т. 2. – М.: ИВМ РАН, 2000. – С. 160 – 188. 35. Левченко О.Ю. Математическое моделирование противобактериального иммунного ответа // Научный журнал КубГАУ. – 2011. – № 66(02). – URL: http://ej.kubagro.ru/2011/02/pdf/05.pdf (дата обращения: 07.12.15). 36. Летов А.М. Математическая теория процессов управления. – М.: Наука, 1981. – 256 с. 123 37. Луговскова Ю.П. Математическое моделирование оптимальных процессов лечения инфекционных заболеваний: автореф. дис… канд. физ.-мат. наук. – Самара: СамГУ, 2009. – 18 с. 38. Марчук Г.И. Математические модели в иммунологии. – М.: Наука, 1980. – 264 с. 39. Марчук Г.И. Математические модели в иммунологии. 2-е изд. – М.: Наука, 1985. – 240 с. 40. Марчук Г.И. Математические модели в иммунологии. Вычислительные методы и эксперименты. – М.: Наука, 1991. – 304 с. 41. Марчук Г.И., Бербенцова Э.П. Острые пневмонии: иммунология, оценка тяжести, клиника, лечение. – М.: Наука, 1989. – 304 с. 42. Марчук Г.И., Бербенцова Э.П. Хронический бронхит: иммунология, оценка тяжести, клиника, лечение. – М.: Успехи физических наук, 1995. – 478 с. 43. Марчук Г.И., Петров Р.В. Вирусное поражение органа и иммунофизиологические реакции защиты (математическая модель) // Проблемы и перспективы современной иммунологии. Методологический анализ. – Новосибирск: Наука, 1988. – С. 100 – 108. 44. Марчук Г.И., Петров Р.В. Математическая модель противовирусного иммунного ответа // Вычислительные процессы и системы. Вып. 1. – М.: Наука, 1983. – C. 5 – 11. 45. Марчук Г.И., Петров Р.В. Математическое моделирование противовирусного иммунного ответа // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 145 – 154. 46. Марчук Г.И., Погожев И.Б., Зуев С.М. Условия подобия процессов в системах взаимодействующих частиц // Докл. РАН. – 1995. – Т. 345, № 5. – С. 605 – 606. 47. Марчук Г.И., Романюха А.А., Бочаров Г.А. Математическое моделирование противовирусного иммунного ответа при вирусном гепатите B // 124 Математические вопросы кибернетики. Вып. 2 / под ред. С.В. Яблонского. – М.: Наука, 1989. – С. 5 – 70. 48. Молер Р. О математике и статистике в иммунологии // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 33 – 45. 49. Нисевич Н.И., Марчук Г.И., Зубикова И.И., Погожев И.Б. Математическое моделирование вирусного гепатита. – М.: Наука, 1981. – 352 с. 50. От иммунологии к демографии: моделирование иммунной истории жизни / А.А. Романюха, С.Г. Руднев, Т.Е. Санникова, Г.И. Марчук // Геронтология in silico: становление новой дисциплины. Математические модели, анализ данных и вычислительные эксперименты [Электронный ресурс]: сборник науч. тр. / под ред. Г.И. Марчука, В.Н. Анисимова, А.А. Романюхи, А.И. Яшина. 2-е изд. (эл.). – М.: БИНОМ. Лаборатория знаний, 2012. – Разд. IV, гл. 13. – С. 396 – 482. 51. Петров Р.В. Иммунология. – М.: Медицина, 1987. – 416 с. 52. Погожев И.Б. Беседы о подобии процессов в живых системах. – М.: Наука, 1999. – 224 с. 53. Погожев И.Б. Некоторые проблемы приложений математических моделей заболеваний в клинической практике // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 175 – 189. 54. Погожев И.Б. Применение математических моделей заболеваний в клинической практике. – М.: Наука, 1988. – 192 с. 55. Поляк Б.Т. Введение в оптимизацию. – М.: Наука, 1983. – 384 с. 56. Понтрягин Л.С. Математическая теория оптимальных процессов. – М.: Наука, 1978. – 392 с. 57. Ройт А., Бростофф Дж., Мейл Д. Иммунология. – М.: Мир, 2000. – 592 с. 58. Романюха А.А. Сопоставление математической модели заболевания и клинических данных // Математическое моделирование в иммунологии 125 и медицине / под ред. Г.И. Марчука. – Новосибирск: Наука, 1982. – С. 27 – 32. 59. Романюха А.А., Бочаров Г.А. Моделирование отдельных процессов иммунного ответа при вирусном гепатите B на основе математической модели противовирусного иммунного ответа // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 155 – 167. 60. Романюха А.А., Каркач А.С. Индивидуально-ориентированная модель динамики инфекционного процесса в неоднородной популяции // Математическое моделирование. – 2003. – Т. 15, № 5. – С. 95 – 105. 61. Романюха А.А., Руднев С.Г. Вариационный принцип в исследовании противоинфекционного иммунитета на примере пневмонии // Математическое моделирование. – 2001. – Т. 13, № 8. – С. 65 – 84. 62. Романюха А.А., Руднев С.Г. Математическое моделирование иммуновоспалительных процессов в лёгких. Поиск оптимальности // Вычислительная математика и математическое моделирование: труды международной конференции / под ред. В.П. Дымникова. Т. 2. – М.: ИВМ РАН, 2000. – С. 212 – 233. 63. Романюха А.А., Руднев С.Г., Зуев С.М. Анализ данных и моделирование инфекционных заболеваний // Современные проблемы вычислительной математики и математического моделирования. Т. 2. Математическое моделирование / под ред. В.П. Дымникова. – М.: Наука, 2005. – С. 352 – 404. 64. Русаков С.В., Чирков М.В. Идентификация параметров и управление в математических моделях иммунного ответа // Российский журнал биомеханики. – 2014. – Т. 18, № 2. – С. 259 – 269. 65. Русаков С.В., Чирков М.В. Идентификация параметров математической модели влияния иммунотерапии на динамику иммунного ответа // Вестник Пермского университета. Математика. Механика. Информатика. – 2014. – Вып.1 (24). – С. 46 – 50. 126 66. Русаков С.В., Чирков М.В. Математическая модель влияния иммунотерапии на динамику иммунного ответа // Проблемы управления. – 2012. – № 6. – С. 45 – 50. 67. Русаков С.В., Чирков М.В. Управление с обратной связью в классической системе типа «хищник-жертва» // Российский журнал биомеханики. – 2015. – Т. 19, № 1. – С. 65 – 72. 68. Русаков С.В., Чирков М.В. Усиление иммунного ответа при хронической форме заболевания в условиях с неполной информацией // Врачаспирант. – 2015. – № 3.2(70). – С. 269 – 275. 69. Сапин М.Р., Этинген Л.Е. Иммунная система человека. – М.: Медицина, 1996. – 304 с. 70. Саридис Дж. Самоорганизующиеся стохастические системы управления. – М.: Наука, 1980. – 400 с. 71. Сейдж Э.П., Мелса Д.Л. Идентификация систем управления. – М.: Наука, 1974. – 248 с. 72. Сейдж Э.П., Уайт Ч.С. Оптимальное управление системами. – М.: Радио и связь, 1982. – 392 с. 73. Смирнова О.А. Математическое моделирование влияния ионизирующей радиации на иммунную систему млекопитающих // Физика элементарных частиц и атомного ядра. – 1996. – Т. 27. Вып. 1. – С. 243 – 292. 74. Смирнова О.А. Радиация и организм млекопитающих: модельный подход. – М., Ижевск: Регулярная и хаотическая динамика. Институт компьютерных исслед., 2006. – 224 с. 75. Смирнова О.А, Степанова Н.В. Электронное моделирование динамики иммунной реакции // Вестник МГУ. Сер. Физика, астрон. – 1971. – № 5. – С. 520 – 526. 76. Современные методы идентификации систем / под ред. П. Эйкхоффа. – М.: Мир, 1983. – 400 с. 77. Степаненко Р.Н., Скалько Ю.И. Математическая модель усиления иммунного ответа стимулятором антителопродукции // Математические 127 модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 136 – 144. 78. Султонова Ш.Х., Меркулова Н.Н. Исследование математической модели усиления иммунного ответа с применением одношаговых и многошаговых методов // Вестн. Том. гос. ун-та. Математика и механика. – 2013. – № 3(23). – С. 51 – 59. 79. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежёсткие задачи. – М.: Мир, 1990. – 512 с. 80. Ходаковский Н.И., Кузьменко Б.В. Построение модели состояния здорового человека на основе работы иммунной системы // Управляющие системы и машины. – 2014. – № 5. – С. 23 – 28. 81. Холл Дж., Уатт Дж. Современные численные методы решения обыкновенных дифференциальных уравнений. – М.: Мир, 1979. – 312 с. 82. Хофман Дж., Купер-Виллис Э. Симметрия, сложность и устойчивость в сетевой теории иммунной системы // Математические модели в иммунологии и медицине / под ред. Г.И. Марчука, Л.Н. Белых. – М.: Мир, 1986. – С. 110 – 122. 83. Хэйл Дж. Теория функционально-дифференциальных уравнений. – М.: Мир, 1984. – 422 с. 84. Цыпкин Я.З. Основы информационной теории идентификации. – М.: Наука, 1984. – 320 с. 85. Чирков М.В. Идентификация параметров и дискретное управление в моделях инфекционных заболеваний // Вестник молодых учёных ПГНИУ [Электронный ресурс]: сб. науч. тр. / отв. редактор В.А. Бячкова. – Пермь: Перм. гос. нац. исслед. ун-т., 2014. – Вып. 4. – С. 332 – 339. 86. Чирков М.В. Моделирование лечения хронической формы инфекционного заболевания в условиях с неполной информацией // Вестник молодых учёных ПГНИУ [Электронный ресурс]: сб. науч. тр. / отв. редактор В.А. Бячкова. – Пермь: Перм. гос. нац. исслед. ун-т., 2015. – Вып. 5. – С. 164 – 170. 128 87. Чирков М.В. Управление в условиях неопределённости на основе математической модели противовирусного иммунного ответа // Научнотехнический вестник Поволжья. – 2017. № 4. – С. 226 – 229. 88. Эйкхофф П. Основы идентификации систем управления: оценивание параметров и состояния. – М.: Мир, 1975. – 684 с. 89. Эльсгольц Л.Э., Норкин С.Б. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. – М.: Наука, 1971. – 296 с. 90. Bajpai P., Chaturvedi A., Dwivedi A.P. Optimal therapeutic control modeling for immune system response // International Journal of Computer Applications. – 2011. – V. 21, N. 4. – P. 27 – 30. 91. Bell G. Mathematical model of clonal selection and antibody production I // J. Theor. Biol. – 1970. – V. 29. – P. 191 – 232. 92. Bell G. Mathematical model of clonal selection and antibody production II // J. Theor. Biol. – 1972. – V. 33. – P. 339 – 378. 93. Bell G. Prey-predator equations simulating an immune response // Math. Biosci. – 1973. – V. 16. – P. 291 – 314. 94. Bellen A.J., Zennaro M. Numerical solution of delay differential equations to an implicit Runge-Kutta method // Numer. Math. – 1985. – V. 47. – P. 301 – 316. 95. Bocharov G.A., Marchuk G.I., Romanyukha A.A. Numerical solution by LMMs of stiff delay differential systems modeling an immune response // Numer. Math. – 1996. – V. 73. – P. 131 – 148. 96. Bocharov G.A., Romanyukha A.A. Numerical treatment of the parameter identification problem for delay-differential system arising in immune response modeling // Applied Numerical Mathematics. – 1994. – V. 15. – P. 307 – 326. 97. Bodnar M., Forys U. A model of immune system with time-dependent immune reactivity // Nonlinear Analysis: Theory, Methods & Applications. – 2009. – V. 70, N 2. – P. 1049 – 1058. 129 98. Bodnar M., Forys U. Behavior of solutions to Marchuk’s model depending on a time delay // Internat. J. Appl. Math. Comput. Sci. – 2000. – V. 10, N 1. – P. 101 – 116. 99. Bodnar M., Forys U. Periodic dynamics in a model of immune system // Appl. Math. – 2000. – V. 27, N 1. – P. 113 – 126. 100. Bruni C., Giovenco M.A., Koch G., Strom R. A dynamical model of humoral immune response // Math. Biosci. – 1975. – V. 27. – P. 191 – 211. 101. Boer de R.J., Hogeweg P. Memory but no suppression in low-dimensional symmetric idiotypic networks // Bull. Math. Biol. – 1989. – V. 51. – P. 223 – 246. 102. Boer de R.J., Hogeweg P. Stability of symmetric idiotypic networks a critique of Hoffmann’s analysis // Bull. Math. Biol. – 1989. – V. 51. – P. 217 – 222. 103. Butcher J.C. Numerical methods for ordinary differential equations. – Chichester: John Wiley & Sons Ltd, 2008. – 490 p. 104. Cooper-Willis A., Hoffman G.V. Symmetry of effector function in the immune system network // Molecular immunology. – 1983. – V. 20. – P. 865 – 870. 105. Cryer C., Tavernini L. The numerical solution of Volterra functional differential equations by Euler’s method // SIAM J. Numer. Anal. – 1972. – V. 9. – P. 105 – 129. 106. Delisi C. Some mathematical problems on the initiation and regulation of the immune response // Math. Biosci. – 1977. – V. 35. – P. 1 – 26. 107. Dibrov B.F., Lifshits M.A., Volkenstein M.V. Mathematical model of immune process // J. Theor. Biol. – 1977. – V. 65. – P. 609 – 631. 108. Fong T., Di Bisceglie A.M., Biswas R., et. al. High levels of viral replication during acute hepatitis B infection predict progression to chronicity // J. Med. Virol. – 1994. – V. 43. – P. 155 – 158. 109. Forys U. Global analysis of Marchuk’s model in a case of strong immune system // J. Biol. Sys. – 2000. – V. 8, N 4. – P. 331 – 346. 130 110. Forys U. Global analysis of Marchuk’s model in a case of weak immune system // Math. Comp. Model. – 1995. – V. 25. – P. 97 – 106. 111. Forys U. Global analysis of the initial value problem for a system of ODE modeling the immune system after vaccination // Math. Comp. Model. – 1999. – V. 29. – P. 79 – 85. 112. Forys U. Hopf bifurcation in Marchuk’s model of immune reactions // Math. Comp. Model. – 2001. – V. 34. – P. 725 – 735. 113. Forys U. Marchuk’s model of immune system dynamics with application to tumour growth // J. Theor. Medicine. – 2002. – V. 4, N. 1. – P. 85 – 93. 114. Forys U. Mathematical model of an immune system with random time of reaction // Appl. Math. – 1991. – V. 21, N.4. – P. 521 – 536. 115. Forys U., Zolek N. A model of immune system after vaccination // ARI. – 1998. – V. 50, N. 3. – P. 180 – 184. 116. Forys U., Zolek N. Complementary analysis of the initial value problem for a systems of ODE modelling immune system after vaccination // Appl. Math. – 2000. – V. 27, N.1. – P. 103 – 111. 117. Galenbuhr V., Bersini H., Stewart J., Varela F.J. Natural tolerance in a simple immune network // J. Theor. Biol. – 1995. – V. 177. – P. 199 – 213. 118. Gollmann L., Maurer H. Theory and applications of optimal control problems with multiple time-delays // J. Ind. Manag. Optim. – 2014. – V. 10, N. 2. – P. 413 – 441. 119. Gunter N., Hoffman G.W. Qualitative dynamics of network model for regulation of the immune system: a rationale for the IgM to IgG switch // J. Theor. Biol. – 1982. – V. 94. – P. 815 – 855. 120. Hege G.S., Cole G. A mathematical model relating circulating antibody and antibody forming cells // J. Immunol. – 1966. – V. 97. – P. 34 – 40. 121. Hoffman G.W. A theory of regulation and self-nonself discrimination in an immune network // Eur. J. Immunol. – 1975. – V. 5. – P. 638 – 647. 131 122. Hoffman G.W. The application of stability criteria in evaluating network regulation model // Regulation of immune response dynamics. – Boca Raton: CRC Press, 1982. – P. 137 – 162. 123. Jilek M. On contact of immunocompetent cells with antigen (Note on probability model) // Folia Microbiol. – 1971. – V. 16. – P. 83 – 87. 124. Jilek M. The number of immunologically activated cells after repeated immunization (A mathematical model) // Folia Microbiol. – 1971. – V. 16. – P. 12 – 23. 125. Jilek M., Sterzl J. A model of differentiation of immunological component cells // Developmental aspects of antibody formation and structure. – Prague: Academia, 1970. – P. 960 – 981. 126. Marchuk G.I. Mathematical modeling of immune response in infectious diseases. – Dordrecht: Springer Science & Business Media, 2013. – 350 p. 127. Marchuk G.I., Romanyukha A.A., Bocharov G.A. Mathematical model of antiviral immune response. II. Parameter identification for acute viral Hepatitis B // J. Theor. Biol. – 1991. – V. 151, N. 1. – P. 41 – 70. 128. Marchuk G.I., Pogozhev I.B., Zuev S.M. Similarity conditions of the processes in system of interacting particles // Russ. J. Numer. Anal. Math. Modeling. – 1996. – V. 11, N 1. – P. 41 – 47. 129. Mohler R. Bilinear control processes. – N.Y.: Academic Press, 1973. – 224 p. 130. Nesterenko V.G. A concept of immune regulation of somatic cell differention // J. Theor. Biol. – 1984. – V. 107. – P. 443 – 456. 131. Nesterenko V.G. Network interaction and the regulation of the immune response // Folia Biologica. – 1984. – V. 30. – P. 231 – 250. 132. Nesterenko V.G. Role of asymmetry in the immune network // Folia Biologica. – 1986. – V. 32. – P. 256 – 272. 133. Oberle H.J., Pesch H.J. Numerical treatment of delay differential equations by hermite interpolation // Numer. Math. – 1981. – V. 37. – P 235 – 255. 134. Perelson A.S. Modelling viral and immune system dynamics // Nat. Rev. Immunol. – 2002. – V. 2, N. 1. – P. 28 – 36. 132 135. Perelson A.S., Mirmirani M., Oster G.F. Optimal strategies in immunology: I. B-cell differentiation and proliferation // J. Math. Biol. – 1976. – V. 3. – P. 325 – 367. 136. Perelson A.S., Mirmirani M., Oster G.F. Optimal strategies in immunology: II. B-memory cell production // J. Math. Biol. – 1978. – V. 5. – P. 213 – 256. 137. Pimbley G.H. Periodic solutions of predator-prey equations simulating an immune response, I // Math. Biosci. – 1974. – N 20. – P. 27 – 51. 138. Pimbley G.H. Periodic solutions of predator-prey equations simulating an immune response, II // Math. Biosci. – 1974. – N 21. – P. 251 – 277. 139. Richter P.H. A network theory of immune system // Eur. J. Immunol. – 1975. – V. 5. – P. 350 – 354. 140. Stengel R.F., Ghigliazza R.M. Stochastic optimal therapy for enhanced immune response // Math. Biosci. – 2004. – V. 191. – P. 123 – 142. 141. Stengel R.F., Ghigliazza R.M., Kulkarni N.V. Optimal control of innate immune response // Optim. Control Appl. & Meth. – 2002. – V. 23. – P. 91 – 104. 142. Stengel R.F., Ghigliazza R.M., Kulkarni N.V. Optimal enhancement of immune response // Bioinformatics. – 2002. – V. 18, N 9. – P. 1227 – 1235. 143. Tavernini L. One-step method for the numerical solution of Volterra functional differential equations // SIAM J. Numer. Anal. – 1971. – V. 8. – P. 786 – 795. 144. Volkov I.K., Zuev S.M. Conditions for identifiability of autonomous systems of ordinary differential equations // Russ. J. Numer. Anal. Math. Modelling. – 1999. – V. 14, N 6. – P. 513 – 525. 133 Приложение