Оптимальные возмущения систем с запаздывающим аргументом для управления динамикой инфекционных заболеваний на основе многокомпонентных воздействий

Полный текст

1. ВВЕДЕНИЕ Новые направления прикладной математики в различных областях естествознания формируются на основе новых понятий и постулатов, конструктивных принципов и методов исследований. Становление математического моделирования в науках о жизни связано с пионерскими работами В. Вольтерра, А. Н. Колмогорова, Н. Винера и др., в которых были предложены такие основы. В иммунологии, одной из бурно развивающихся дисциплин современной биомедицины, основы математического моделирования были заложены Г. И. Марчуком [2, 4-6, 33, 34]. Им был предложен и успешно реализован подход, связанный с «координатизацией» (понятие, введенное Г. Вейлем [40]) явлений сложной природы, изучаемых иммунологами, таких как тяжесть заболевания, уровень иммунитета, иммунореактивность. В основу построения проблемно-ориентированных математических моделей иммунных процессов был положен принцип выделения из большого числа признаков и факторов изучаемых процессов набора «главных» параметров, доступных количественному оцениванию по клиническим и экспериментальным данным. Был найден некоторый «оптимальный» уровень идеализации для математического описания сложной системы на основе нелинейных дифференциальных уравнений с запаздывающим аргументом, в рамках которого оказалось возможным проводить калибровку и исследование моделей для реальных инфекций, а также строить логические выводы, имеющие содержательную ценность для иммунологии. Г. И. Марчук задолго Qc РОССИЙСКИЙ УНИВЕРСИТЕТ ДРУЖБЫ НАРОДОВ, 2017 392 до появления современной системно-биологической парадигмы исследования иммунных процессов [19, 27, 29] осознал необходимость применения многопараметрических моделей иммунного ответа: «... будущее медицины - лечение индивидуального больного на основе слежения за его индивидуальными иммунологическими, эндокринологическими, сосудистыми особенностями с учетом непрерывно приобретаемых с возрастом хронических локусов различной этиологии...» на основе многомасштабных моделей. Такой взгляд полностью отвечает современному развитию системной иммунологии [20, 31], которое нацелено на разработку и применение методов коррекции иммунофизиологических процессов на различных уровнях регуляции (генно-молекулярном, клеточном, тканевом, органном и системном). Таким образом, изучение многоуровневых регуляторных механизмов, лежащих в основе наблюдаемого спектра динамики вирусных заболеваний в категориях причинно-следственных взаимоотношений и откликов организма человека на действие противовирусных и иммуномодулирующих препаратов, является фундаментальной междисциплинарной задачей. Разработка математических моделей для описания и исследования динамики инфекционных заболеваний имеет своей конечной целью анализ характеристик чувствительности процессов к воздействиям различной природы, изменяющих параметры соответствующей динамической системы и ее состояние. На их основе могут быть рассчитаны количественные характеристики устойчивости и чувствительности фенотипически различных траекторий динамики системы. В рамках реализуемого в настоящее время системного подхода мы рассматриваем инфекционное заболевание как траекторию некоторой динамической системы в пространстве признаков. Неблагоприятное течение заболевания можно рассматривать как выход за пределы области в окрестности «нормальной» траектории или аттрактора. Математическая модель инфекционного заболевания используется для исследования характеристик робастности и хрупкости динамики системы [28] по отношению к возмущениям параметров процессов или состояния системы с целью поиска набора управляющих воздействий для коррекции неблагоприятной динамики. В работе [12] впервые для задач математической иммунологии нами был реализован подход к построению возмущений начального состояния динамической системы с запаздыванием с помощью так называемых оптимальных возмущений, обеспечивающих максимально возможную амплификацию возмущения в заданной норме. В основе этого подхода - методы, заимствованные из теории аэродинамической устойчивости и обобщенные на системы с запаздыванием. С помощью оптимальных возмущений, выбираемых из подпространства кусочно-постоянных функций, были найдены малые возмущения стационарного состояния системы с запаздыванием, оказывающие максимальное воздействие на ее динамику с точки зрения заданных критериев. Для двух типов устойчивых стационарных состояний модели экспериментальной вирусной инфекции, отвечающих биологическим сценариям персистенции вирусов ниже порога обнаружения и хроническим вирусным инфекциям с высоким уровнем вирусной нагрузки, были рассчитаны несколько структурно различных оптимальных возмущений и исследовано их влияние на динамику инфекционного процесса. При этом были сформулированы следующие вопросы, требующие дальнейшего исследования, которым и посвящена данная работа: 1. Значимость отдельных компонент оптимальных возмущений или, в более широком плане, построение наиболее простых оптимальных возмущений при заданных ограничениях на длительность возмущения и эффективность воздействия. 2. Выбор подпространства, в котором ищется возмущение, и его связь с характеристиками реакции инфекционного и иммунного процессов на действие антивирусных и иммуномодулирующих препаратов с учетом их фармакокинетики и фармакодинамики. Работа состоит из семи разделов. В разделе 2 кратко поясняется, что такое оптимальные возмущения динамических систем без запаздывания и обсуждаются известные методы их вычисления. В разделе 3 описана математическая модель динамики экспериментальной вирусной инфекции, сформулированная в виде системы четырех нелинейных дифференциальных уравнений с запаздывающим аргументом. В модели учитываются два различных запаздывания, отражающие длительность деления клеток иммунной системы в ходе формирования противовирусного иммунного ответа и длительность программирования процесса гибели клеток, вызванного присутствием вирусных антигенов. Исследуются два биологически содержательных положения равновесия системы, с низкой и высокой вирусной нагрузкой. Они соответствуют различным вариантам хронического вирусного заболевания. Рассматриваются общие закономерности фармакокинетики и фармакодинамики противовирусных и иммуномодулирующих препаратов, которые определяют выбор подпространства начальных возмущений для задачи управления. В разделе 4, следуя работе [12], вводится понятие оптимального возмущения для системы с запаздыванием, а в разделе 5 описывается простейший метод его вычисления. Результаты расчетов оптимальных возмущений при различных значениях их параметров приводятся в разделе 6, где также дается их содержательная интерпретация. Итог работы подводится в заключительном разделе 7. 1. ОПТИМАЛЬНЫЕ ВОЗМУЩЕНИЯ СТАЦИОНАРНЫХ РЕШЕНИЙ ДИНАМИЧЕСКИХ СИСТЕМ В качестве возмущений стационарного состояния системы мы будем использовать так называемые оптимальные возмущения, широко применяемые в теории аэродинамической устойчивости в рамках моделей без запаздывания. В аэродинамике различают естественный (надкритический) и обходной (докритический) сценарии ламинарно-турбулентного перехода [13]. С увеличением числа Рейнольдса при превышении им критического значения, определяющего границу устойчивости к бесконечно малым возмущениям, ламинарное течение заведомо теряет устойчивость, что, как правило, приводит к его турбулизации. При этом основную роль играют неустойчивые моды. Однако из-за наличия в реальных течениях малых конечных возмущений ламинарно-турбулентный переход на практике часто происходит при докритических числах Рейнольдса (докритические сценарии перехода). Одним из основных факторов, вызывающих докритический ламинарно-турбулентный переход, является возможность существенного роста энергии возмущения на конечных временных интервалах. Этот рост обеспечивают возмущения, представляющие собой суперпозиции большого числа существенно неортогональных устойчивых мод. Развитие такого малого возмущения может привести к переходу основного течения в квазистационарное линейно неустойчивое состояние, в котором начинает развиваться вторичная неустойчивость, приводящая к ламинарно-турбулентному переходу. Максимально возможный при данном докритическом числе Рейнольдса «подскок» (v(t)) Γmax = sup E E(v(0)) средней плотности кинетической энергии возмущения, где v(t) означает скорость возмущения в момент времени t и супремум берется по всем допустимым ненулевым начальным возмущениям и всем t > 0, называют максимальной амплификацией средней плотности кинетической энергии возмущений. Начальное возмущение v(0), при котором достигается Γmax, называют оптимальным возмущением. Задача определения максимальной амплификации средней плотности кинетической энергии возмущений в рамках полных уравнений эволюции возмущений представляет значительную сложность, и в настоящее время нет каких-либо известных алгоритмов ее решения. Однако для моделирования и анализа докритических сценариев ламинарно-турбулентного перехода достаточно рассматривать возмущения, оптимальные в рамках линеаризованных уравнений [24, 38]. Поэтому на практике ограничиваются вычислением максимальной амплификации и соответствующих оптимальных возмущений для линеаризованных уравнений. Поскольку для линеаризованных уравнений максимальная амплификация не зависит от величины первоначального возмущения, для определенности полагают E(v(0)) = 1. Оптимальные возмущения могут играть существенную роль и в поведении любой другой динамической системы в окрестности положения равновесия, если уравнения эволюции возмущений этой системы, линеаризованные относительно стационарного состояния, имеют существенно неортогональные моды, отвечающие существенно различным собственным значениям (этому условию удовлетворяют, например, модели ядерных реакторов с запаздывающими нейтронами [3, 30]). В этом случае всегда найдутся начальные возмущения, нормы которых будут значительно возрастать на конечных временных интервалах. Такие возмущения могут существенно повлиять на динамику системы. Следует отметить, что не только оптимальные возмущения, но вообще возмущения, растущие по норме на конечных временных интервалах, весьма специфичны. Как правило, если начальное возмущение устойчивого стационарного состояния выбирать случайным образом, то норма возмущения будет монотонно убывать с вероятностью единица. Вычисление оптимального возмущения требует специальных алгоритмов. Разработанные ранее методы [14, 15, 36, 37] позволяют эффективно вычислять максимальную амплификацию второй (евклидовой) нормы решения и оптимальные возмущения с заданной точностью для систем обыкновенных дифференциальных уравнений без запаздывания. Пусть v означает возмущение устойчивого стационарного состояния некоторой конечномерной автономной динамической системы и линеаризованные уравнения эволюции возмущения этого стационарного состояния записаны в виде dv = Av, (2.1) dt где A - заданная постоянная матрица. Тогда v(t) = exp{tA}v(0) и расчет максимальной амплификации второй нормы возмущений сводится к вычислению максимальной амплификации второй нормы решений задач Коши для системы (2.1), которая очевидно равна t�0 Γmax = max exp{tA} 2. (2.2) При этом оптимальным возмущением будет правый сингулярный вектор vopt матрицы exp{toptA}, отвечающий ее максимальному сингулярному числу [21], где topt означает значение t, при котором достигается максимум (2.2). Обычно множество оптимальных возмущений ограничивают теми из них, при которых максимальная амплификация достигается при наименьшем значении t, то есть при t�0 topt = min arg max exp{tA} 2. Описанный выше подход к вычислению максимальной амплификации второй нормы возмущений и соответствующих оптимальных возмущений сводится, таким образом, к задаче вычисления для заданной квадратной комплексной матрицы A глобального максимума функции Γ(t) = exp{tA} 2, и минимального topt, при котором достигается этот максимум. Предполагается, что спектр матрицы A лежит в левой полуплоскости (условие устойчивости рассматриваемого стационарного состояния), а максимальное собственное значение hmax матрицы A + A∗ положительное (иначе topt = 0, а Γmax = 1). Найти tˆopt, дающее максимум Γ(t) с заданной относительной точностью δ, можно, вычислив Γ(t) на равномерной сетке с узлами tj = τj и достаточно мелким шагом τ. Для этого вычисляем матрицу E1 = exp{τ A} методом, описанным в [25]. В остальных узлах сетки значение матричной экспоненты Ej = exp{tj A} вычисляем по формуле Ej = E1Ej-1. Нормы полученных матриц вычисляем с помощью соответствующей процедуры пакета LAPACK [9]. Остановить вычисления достаточно, если Ej 2 < 1 (можно показать, что в этом случае tj > topt), взяв в качестве tˆopt узел сетки, в котором норма матричной экспоненты максимальная. Этот алгоритм мы далее будем называть простейшим. Основываясь на неравенстве Γ(t)2 � exp(hmaxt), непосредственно следующем из очевидной оценки d 2 2 dt v(t) 2 � hmax v(t) 2, несложно показать, что найденное таким образом tˆopt будет удовлетворять неравенству |Γ(tˆopt)2 - Γ(topt)2|/Γ(topt)2 � δ, (2.3) если выбрать τ = ln(1 + δ)/hmax. Однако, как правило, вычислительные затраты такого алгоритма будут катастрофически большими, поскольку из-за большой величины topthmax будет очень большой величина j, при которой Γ(tj ) достигает максимума. Для уменьшения вычислительных затрат в работе [36] было предложено заменять матрицу A ее формой Шура [21] (то есть верхней треугольной матрицей, унитарно подобной исходной) с диагональными элементами, упорядоченными по невозрастанию вещественных частей, а для ускорения вычисления матриц Ej, которые для верхней треугольной матрицы будут также верхними треугольными, и их норм применять аппроксимацию вида Ek ≈ Pk 0 Qk, (2.4) где Pk - квадратная верхняя треугольная матрица меньшего порядка, чем Ek, а Qk - унитарная прямоугольная матрица. Учитывая, что норма правой части (2.4) равна Pk 2, такая аппроксиj-k мация позволит вычислять, начиная с j = k + 1, матрицы E�1 Pk вместо Ej, где E�1 означает главную подматрицу матрицы E1 того же порядка, что и Pk. В работе [36] показано, что если такую редукцию выполнять время от времени в процессе работы исходного алгоритма, то вычислительные затраты существенно уменьшатся (в тысячи раз) при той же точности результата в смысле (2.3). Аппроксимацию (2.4) можно выполнить с заданной точностью следующим образом. Разбиваем матрицу Ek по строкам на два блока: Ek = E Г (1) l k E (2) k так, чтобы норма нижнего блока не превосходила заданную малую пороговую величину, и, отбросив нижний блок, выполняем разложение верхнего блока в произведение квадратной верхней треугольной и унитарной прямоугольной матриц, что сводится к вычислению QL-разложения транспонированной матрицы. В пакете LAPACK имеются процедуры, необходимые для вычисления всех упомянутых выше разложений. Поэтому вычисление максимальной амплификации средней плотности кинетической энергии возмущений и соответствующего оптимального возмущения можно организовать так, чтобы основной объем вычислений выполнялся стандартными матричными процедурами. Этот, а также методы, предложенные в недавней работе [37] для систем с большими разреженными матрицами и основанные на минимизации функционала (exp{tA}v0, exp{tA}v0) (v0, v0) по t и v0, представляется целесообразным в будущем перенести на системы с запаздыванием. В данной работе мы использовали алгоритм, предложенный в работе [12] и являющийся обобщением на системы с запаздыванием описанного выше простейшего алгоритма. 2. МАТЕМАТИЧЕСКИЕ МОДЕЛИ ВИРУСНЫХ ИНФЕКЦИЙ И ФАРМАКОКИНЕТИКИ ЛЕЧЕБНЫХ ПРЕПАРАТОВ Математические модели вирусных инфекций описывают процессы деления, дифференцировки и гибели клеток иммунной системы, которые характеризуются значительной продолжительностью по сравнению с характерными временами размножения вирусов. По этой причине использование дифференциальных уравнений с запаздывающим аргументом является адекватным способом моделирования иммунных процессов [2, 4, 5, 33]. Модели описывают различные варианты динамики вирусных инфекций, включая стационарные решения. Стационарные решения соответствуют хроническим вариантам течения заболеваний и являются объектом изучения для разработки новых подходов к их лечению. В основе таких подходов лежит комбинация противовирусных и иммуномодулирующих препаратов. Для теоретического исследования возможностей построения эффективных мультимодальных воздействий на хронические вирусные инфекции нами рассмотрена математическая модель динамики инфекции, вызванной вирусами лимфоцитарного хориоменингита (ВЛХМ), разработанная в [10]. Данная модель сформулирована в виде системы нелинейных дифференциальных уравнений с запаздывающим аргументом и описывает динамику следующих переменных, зависящих от времени: концентрации вирусных частиц V (t), двух популяций специфичных к ВЛХМ цитотоксических лимфоцитов (ЦТЛ) - клеток-прекурсоров Ep(t) и клетокэффекторов Ee(t), а также кумулятивной вирусной нагрузки W (t). Система уравнений имеет вид: d - V (t) = βV (t) 1 dt d 0 V (t) Vmvc § γVEEe(t)V (t), dt Ep(t) = αEp (Ep - Ep(t)) + βpgp(W )V (t - τ )Ep(t - τ ) - - αAP V (t - τA)V (t)Ep(t), (3.1) d dt Ee(t) = bdge(W )V (t - τ )Ep(t - τ ) - - αAEV (t - τA)V (t)Ee(t) - αEe Ee(t), d dt W (t) = bW V (t) - αW W (t), где gp(W ) = 1/(1 + W/θp)2, ge(W ) = 1/(1 + W/θE )2. Биологический смысл параметров модели пояснен в табл. 1. ТАБЛИЦА 1. Биологический смысл параметров модели (3.1). Параметр Биологические смысл β γVE Vmvc τ bp bd θp θE αEp αEe E0 p τA αAP αAE bW αW Константа скорости репликации вирусных частиц Константа скорости элиминации вирусов за счет клеток-эффекторов Максимально возможная концентрация вирусных частиц в селезенке Характерная продолжительность цикла деления ЦТЛ Константа скорости стимуляции ЦТЛ Константа скорости дифференциации ЦТЛ Порог вирусной нагрузки для перехода прекурсоров в состояние анергии Порог вирусной нагрузки для перехода эффекторов в состояние анергии Константа скорости естественной гибели прекурсоров Константа скорости естественной гибели эффекторов Концентрация прекурсоров в селезенке мыши, не имевшей контакта с ВЛХМ Характерная продолжительность перехода ЦТЛ к апоптозу Константа скорости апоптоза прекурсоров Константа скорости апоптоза эффекторов Константа скорости роста кумулятивной вирусной нагрузки Константа скорости восстановления организма от воздействий вирусной нагрузки Для определения решения при t > 0 достаточно задать значения V (t) при -τA � t � 0, значения Ep(t) при -τ � t � 0, значения Ee(0) и W (0). Однако для единообразия далее мы будем предполагать, что начальные значения для всех переменных модели заданы при -τA � t � 0. Задача Коши для системы уравнений (3.1) с неотрицательными начальными условиями и неотрицательными параметрами глобально разрешима на любом конечном интервале времени t ∈ [0,T ]. Данное утверждение можно доказать, используя технику, описанную в [4], на основе метода шагов Беллмана [1], с использованием теоремы Такуэти о положительной инвариантности систем обыкновенных дифференциальных уравнений (ОДУ) [39] и рассматривая линейную систему ОДУ, мажорирующую правую часть системы уравнений модели. Обозначив вектор переменных системы (3.1) через U (t) = (V (t), Ep(t), Ee(t),W (t))T , (3.2) будем записывать эту систему в следующем компактном виде: d dt U (t) = F (U (t),U (t - τ ),U (t - τA)). (3.3) В соответствии со сказанным выше будем считать, что вектор переменных модели U (t) задан при -τA � t � 0. Модель (3.3) допускает существование различных стационарных решений, которым можно дать различную биологическую интерпретацию. Нами рассмотрены два варианта стационарных решений модели. Первое стационарное состояние соответствует латентной форме инфекции с низким уровнем вирусной нагрузки и большой численностью клеток иммунологической памяти [17]. Второе стационарное решение описывает клинически выраженную хроническую инфекцию с большой вирусной нагрузкой и частичным истощением популяции Т-лимфоцитов, специфичных к антигенам вируса [18, 35]. Вариацией параметров в области допустимых значений, указанных в [32], были найдены стационарные состояния, описывающие два вышеуказанных варианта хронической инфекции. Стационарные состояния вычислялись на основе системы нелинейных алгебраических уравнений модели F (U, U, U ) = 0 с помощью метода Ньютона. Поиск наборов параметров, при которых стационарные состояния модели (3.3) обладают требуемыми свойствами, и соответствующих начальных значений для метода Ньютона был выполнен, основываясь на результатах численного бифуркационного анализа, проведенного в работе [32]. Были определены два набора параметров, для каждого из которых было найдено устойчивое положение равновесия. Значения параметров модели (3.3), соответствующие указанным двум стационарным состояниям UI и UII, приведены в табл. 2, а значения компонент этих стационарных состояний - в табл. 3. Соответственно, целью возмущения системы в первом случае может являться либо активация инфекционного процесса с последующим удалением вирусного резервуара, либо непосредственная нейтрализация инфекции, а во втором - восстановление реактивности истощенного звена иммунной системы и снижение вирусной нагрузки. Заметим, что данные сценарии динамики процесса имеют место при ВИЧ инфекции и соответствуют фенотипически различным вариантам инфекции, наблюдаемым у «элит-контроллеров» и «прогрессоров» соответственно [22, 26]. Для лечения хронических форм инфекционных заболеваний перспективным представляется применение комбинированной терапии, т. е. сочетание антиретровирусных и иммуномодулирующих воздействий. Результатом подобных возмущений системы вирус-организм хозяина будет изменение динамики системы. Нас интересует построение мультикомпонентных воздействий на состояние системы в положении равновесия, которые будут иметь своим следствием формирование реакции, характеризующейся активацией иммунного ответа и уменьшением численности вирусной популяции. Для параметризации эффекта комбинированной терапии воспользуемся базовыми зависимостями фармакокинетики, описывающими распространение препаратов в организме. Если в момент времени t0 была введена некоторая доза препарата, то концентрация препарата C(t, t0) тождественно равна нулю при t < t0 и меняется далее по закону, который определяется соответствующей фармакокинетической моделью [8], в том числе, C(t, t0) = a1 exp{-α1(t - t0)} в случае однокамерной модели, a1 exp{-α1(t - t0)}- a2 exp{-α2(t - t0)} в случае однокамерной модели с всасыванием, a1 exp{-α1(t - t0)} + a3 exp{-α3(t - t0)} в случае двухкамерной модели и a1 exp{-α1(t - t0)}- a2 exp{-α2(t - t0)} + a3 exp{-α3(t - t0)} в случае двухкамерной модели с всасыванием, где aj, αj означают некоторые положительные величины. Если за время лечения пациент принимал лекарства l раз в моменты времени t1,..., tl , то c учетом конкретного выбора 0 0 терапии концентрация препарата в любой момент времени будет равна l 0 C(t) = \ Ci(t, ti ), i=1 где Ci(t, ti ) означает концентрацию препарата, введенного в момент времени ti . 0 0 В дополнение к приведенным стандартным моделям фармакокинетики, следует рассматривать и модели, описывающие процессы увеличения численностей клеток и антигена, реализуемых путем адаптивного переноса вирус-специфических лимфоцитов (введение в организм реципиента генетически близких лимфоцитов, которые продолжают функционировать в организме реципиента) и введения вирусных антигенов, соответственно. С учетом поведения перечисленных выше фармаконетических моделей, в качестве базисных функций начального возмущения представляется корректным использовать функции, качественно аппроксимирующие поведение препаратов в рамках однокамерных и двухкамерных моделей, в том числе, и с положительными показателями экспонент, а именно, функции вида (0, 0 � t < t0, для переменных V и W и exp{β(t - t0)}- 1, t0 � t � 0, (0, 0 � t < t0, exp{β(t - t0)}, t0 � t � 0, (3.4) (3.5) для переменных Ep и Ee. Примеры таких функций приведены на рис. 1-2. Заметим, что линейная комбинация нескольких таких воздействий может приводить к немонотонным по времени возмущениям стационарного состояния системы. 3. ОПТИМАЛЬНЫЕ ВОЗМУЩЕНИЯ СТАЦИОНАРНЫХ СОСТОЯНИЙ ДИНАМИЧЕСКИХ СИСТЕМ С ЗАПАЗДЫВАНИЕМ В работе [12] для систем c запаздыванием были впервые предложены физически обоснованные аналоги максимальной амплификации нормы решения и оптимального возмущения. Описанию этих новых понятий посвящен данный раздел. В отличие от работы [12], мы введем более общие определения, позволяющее варьировать значимость отдельных компонент отклика на оптимальное возмущение. Нас будет интересовать поведение системы (3.3) вблизи устойчивого стационарного состояния. Записывая решение вблизи стационарного состояния в виде U (t) = U + εU ∗(t) + O(ε2), где ε - малый по абсолютной величине параметр, подставляя это решение в (3.3) и требуя, чтобы полученные равенства выполнялись при всех сколь угодно малых значениях ε, для функции U ∗(t) получим следующую систему линейных дифференциальных уравнений: d где dt U ∗(t) = L0U ∗(t)+ Lτ U ∗(t - τ )+ LτA U ∗(t - τA), (4.1) ⎛ 2βV ⎞ V β - - γVEEe 0 -γVEV 0 ⎜ vmc ⎟ ⎟ 2 ⎜ ⎜ -2bpV Epgp(W ) ⎟ ⎜ , L0 = ⎜ -αAP V Ep -αEp - αAP V 0 ⎟ ⎟ θp + W ⎟ ⎜ ⎜ 2 -2bdV Epge(W ) ⎟ - ⎜ αAE ⎝ V Ee 0 -αEe § αAEV ⎟ θE + W ⎠ bW 0 0 -αW ⎛ 0 0 0 0⎞ ⎛ 0 0 0 0⎞ Lτ = ⎜bpEpgp(W ) bpVgp(W ) 0 0⎟ , Lτ ⎜-αAP V Ep 0 0 0⎟ . ⎜ ⎟ A = ⎜ ⎟ ⎝bdEpge(W ) bdV ge(W ) 0 0⎠ 0 0 0 0 ⎝-αAEV Ee 0 0 0⎠ 0 0 0 0 Систему (4.1) будем называть линеаризованными уравнениями эволюции возмущений. Начальные значения этой системы будем задавать, как и начальные значения системы (3.3), в интервале -τA � t � 0. Для вектора переменных системы (4.1) введем в рассмотрение следующее семейство локальных норм в момент времени t: r ⎛ t ⎞1/2 2 U ∗ D,t = ⎝ t-τA DU ∗(ξ) 2 dξ⎠ , (4.2) где D - заданная положительно определенная диагональная матрица, а . 2 - вторая (евклидова) векторная норма. opt Под оптимальным возмущением U ∗ (t) будем понимать решение системы (4.1), обеспечивающее максимальный подскок значения локальной нормы (4.2) по сравнению с ее первоначальным значением, то есть такое возмущение U ∗(t), при котором достигается максимум величины max U ∗ SD,t , t�0 U ∗ D,0 где S - заданная положительно определенная диагональная матрица, позволяющая варьировать значимость отдельных компонент отклика на возмущение. Поскольку оптимальное возмущение по определению является некоторым решением линейной системы (4.1) и, следовательно, полностью определяется своим значением в интервале -τA � t � 0, при построении оптимального возмущения наряду с выбором нормы, в которой проводится оптимизация, принципиальным является вопрос о том, из какого пространства функций, заданных в интервале -τA � t � 0, мы берем начальные значения. Это подпространство функций [-τA, 0] → R4 мы далее будем обозначать через Q. Для существования максимума необходимо, чтобы Q было полным пространством в норме . D,0. На практике достаточно выбирать в качестве Q линейную оболочку какого-либо конечного набора базисных функций. Это, в частности, гарантирует его полноту. Искать оптимальное возмущение удобнее в два этапа. Сначала вычисляем максимальную амплификацию Γ(t) = max U ∗ SD,t (4.3) Ut U ∗ D,0 локальной нормы решения системы (4.1), где максимум берется по всем решениям, начальные значения которых ненулевые и принадлежат Q, и находим t = topt, при котором функция Γ(t) достигает максимального значения. Если таких t несколько, то для определенности берем из них минимальное. Таким образом, Затем находим topt = min arg max Γ(t). t�0 ∗ U ∗ Uopt ∈ arg max SD,topt . (4.4) Ut U ∗ D,0 Если D, S и Q фиксированы, то любое найденное оптимальное возмущение обеспечивает один и тот же максимальный подскок локальной нормы решения. Обычно максимальная амплификация имеет только одну точку максимума, а решение оптимизационной задачи (4.4) однозначно с точностью до ненулевой мультипликативной константы. Найденное оптимальное возмущение мы будем использовать для возмущения стационарных состояний исходной нелинейной модели (3.3). Для этого в качестве начальных значений мы будем брать opt U (t) = U + εU� ∗ (t) (4.5) opt при -τA � t � 0, где U� ∗ (t) означает нормированное оптимальное возмущение, а ε - параметр. Варьируя этот параметр, можно увеличивать или уменьшать начальное возмущение стационарного состояния. Если абсолютная величина ε мала, то следует ожидать, что полученное решение U (t) системы (3.3) будет близко к (4.5) при t > 0. При больших по абсолютной величине ε из-за влияния нелинейности решение системы (3.3) будет значительно отличаться от (4.5) при t > 0. Существенную роль играет и знак ε. В зависимости от него плотность популяции вируса при t > 0 либо начинает нарастать, либо убывать. Если оптимальное возмущение в (4.5) нормировано так, что первая компонента вектора opt L0U� ∗ opt (0) + Lτ U� ∗ (-τ )+ Lτ A U� ∗ opt (-τA) является положительной, то плотность популяции вируса нарастает при t > 0 в случае ε > 0. 4. ВЫЧИСЛЕНИЕ ОПТИМАЛЬНЫХ ВОЗМУЩЕНИЙ Оптимальные возмущения можно вычислять на основе любой разностной схемы, подходящей для решения задач Коши для систем линейных обыкновенных дифференциальных уравнений с запаздывающим аргументом. Следуя работе [12], мы будем использовать неявную схему второго порядка BDF2 [23] на равномерной сетке {tk = δk : k = -mA + 1, -mA + 2,.. .}, построенной в полуинтервале (-τA, ∞) с шагом δ > 0. Величины m = [τ /δ] и mA = [τA/δ], где [.] означает целую часть числа, являются дискретными аналогами задержек τ и τA соответственно. После дискретизации описанным выше способом система (4.1) примет вид 1,5Uk - 2Uk-1 + 0,5Uk-2 = L U + L U + L U , k = 1, 2,..., (5.1) δ 0 k τ k-m τA k-mA где Uk - сеточная функция, аппроксимирующая U (tk ). В качестве начальных данных для решения задачи Коши нужно задать значения U-mA+1,..., U0. Запишем уравнение (5.1) в виде Uk = C1Uk-1 + C2Uk-2 + CmUk-m + CmA Uk-mA , (5.2) где C1 = 2(1,5I - δL0)-1, C2 = -0,5(1,5I - δL0)-1, A Cm = (1,5I - δL0)-1δLτ , Cm A = (1,5I - δL0)-1δLτ , а I означает единичную матрицу четвертого порядка, и дополним уравнение (5.2) тождествами Uj = Uj, j = k - 1,...,k - mA + 1. Полученную таким образом систему из mA уравнений можно записать в виде где Xk = MXk-1, (5.3) ⎛ Uk ⎞ ⎛ M11 ··· M1mA ⎞ Xk = ⎜ .. ⎟ ⎜ .. .. ⎟ ⎝ . ⎠ , M = ⎝ . . ⎠ . (5.4) Uk-mA+1 MmA1 ··· MmAmA Матрица M в (5.4) является блочной, блочного порядка mA с блоками порядка 4. Все блоки этой матрицы нулевые, кроме поддиагональных блоков Mj+1,j = I (j = 1,..., mA - 1) и четырех блоков M11 = C1, M12 = C2, M1m = Cm и M1mA = CmA , стоящих в первой блочной строке. В силу (5.3), (5.4) сеточный аналог Γk максимальной амплификации (4.3) нормы решения можно записать следующим образом: Γk = max X0∈span(Q)\{0} , H� Mk X0 2 HX0 2 где Q - матрица размера 4mA × p (p � 4mA), столбцы которой образуют базис в сеточном аналоге подпространства Q, span(Q) означает линейную оболочку столбцов матрицы Q, H = ImA ⊗ D, H� = ImA ⊗ (SD), D и S - диагональные матрицы, задающие локальные нормы, в которых вычисляется оптимальное возмущение, а ⊗ - символ Кронекерова произведения. С учетом того, что Γk = H� Mk H-1Q� 2, где Q� - матрица, полученная ортогонализацией столбцов матрицы HQ, вычисление Γk сводится к формированию матрицы Y0 = H-1Q� и вычислению матриц Yk по рекуррентной формуле Yk = MYk-1 с одновременным вычислением норм матриц H� Yk. Пусть kopt означает значение k, при котором достигается максимум Γk. Вычислив нормированный правый сингулярный вектор η матрицы H� Mkopt H-1Q�, (5.5) 0 отвечающий ее максимальному сингулярному числу [21], начальное значение Xopt сеточного аналога оптимального возмущения Xopt можно найти по формуле Xopt = H-1Q�η. k 0 Отметим, что для повышения эффективности описанного выше алгоритма матрицу M надо хранить и умножать на векторы в разреженном формате. 5. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ Оптимальные возмущения искались в подпространстве Q заданных в интервале [0, τA] функций U (t), представляющих собой экспоненциальные зависимости вида (3.4) для переменных V и W и вида (3.5) для переменных Ee и Ep. В качестве значений параметра t0 для этих базисных функций брались величины -iτA/l (i = 0,...,l - 1), то есть правые концы l равных подынтервалов, на которые предварительно разбивался интервал [-τA, 0]. В качестве весов нормы (диагональные элементы матрицы D) брались величины, обратные к компонентам рассматриваемого стационарного решения. Для расчета оптимальных возмущений использовалась сетка с шагом δ = 5 · 10-3. Для стационарного состояния UI вычислялось одно оптимальное возмущение, соответствующее следующему набору параметров базисных функций: l = 2, β = -1, 0,5 и 1. Матрица S была взята единичной, т. е. все компоненты возмущения считались равнозначными. Далее это возмущение мы будем называть первым. Для стационарного состояния UII были найдены и исследованы четыре оптимальных возмущения, которые мы далее будем называть вторым, третьим, четвертым и пятым, соответственно. При расчете второго оптимального возмущения в качестве Q выбиралось подпространство функций, у которых компонента W является тождественным нулем. При расчете остальных оптимальных возмущений ненулевыми допускались все компоненты. Параметры возмущений выбирались следующими: для второго оптимального возмущения l = 2, β = -1 для компоненты V и l = 14, β = 1 для компонент Ep и Ee; для третьего и четвертого оптимальных возмущений l = 2, β = -1, 0,5 и 1 для всех компонент; для пятого оптимального возмущения l = 2, β = -12, -6, -1, 0,5 и 1 для всех компонент. При расчете второго и четвертого оптимальных возмущений матрица S была выбрана единичной, при расчете третьего и пятого оптимальных возмущений диагональные элементы этой матрицы, отвечающие компонентам Ep и Ee, были увеличены в 10 раз, что повысило значимость этих компонент. Во всех рассмотренных случаях максимальная амплификация, т. е. максимальное сингулярное число Γ(topt) матрицы (5.5), и второе максимальное сингулярное число Γ�(topt) этой матрицы были достаточно хорошо отделены друг от друга (см. табл. 4), что свидетельствовало о единственности оптимального возмущения с точностью до ненулевой мультипликативной константы. Для интегрирования системы (3.3) с начальными данными (4.5) шаг сетки выбирался равным 10-4 и использовалась та же схема BDF2, что и при вычислении оптимальных возмущений. Оптимальные возмущения интерполировались на более мелкую сетку с помощью стандартного алгоритма кубической интерполяции, сохраняющего монотонность и выпуклость. Параметр ε выбирался по абсолютной величине таким, чтобы получить достаточно большой отклик при малом начальном возмущении. Отметим, что для численного интегрирования исходной системы (3.1) ее авторы использовали код DIFSUB-DDE, предназначенный для численного решения жестких систем нелинейных дифференциальных уравнений с постоянными запаздываниями [11]. В коде реализована модификация Гира линейного многошагового метода по схеме BDFp переменного порядка p � 6 с обработкой точек разрыва производных до p +1 порядка, при этом для аппроксимации запаздывающих компонент используется интерполяционный вектор Нордсика. С помощью этого кода были проведены верификационные расчеты, показавшие хорошее совпадение с результатами расчетов с помощью метода BDF2. Стационарное состояние UI, возмущенное начальное значение, полученное добавлением к этому стационарному состоянию оптимального возмущения, взятого с коэффициентом ε, и решения двух задач Коши с такими начальным значением приведены на рис. 3 и 4. Горизонтальные линии соответствуют стационарному состоянию, темные на четырех верхних картинках - возмущенному начальному значению, а темные и светлые на четырех нижних картинках - решениям задач Коши, полученным соответственно с помощью полных и линеаризованных уравнений. В верхней части табл. 3 приведены округленные до трех значащих десятичных цифр максимальные и минимальные значения компонент решений, полученных с помощью полных уравнений. Аналогичные результаты для стационарного состояния UII представлены на рис. 5-8 и в нижней части табл. 3. Для стационарного состояния UI с низкой вирусной нагрузкой целью возмущения являлось обострение инфекционного процесса и выведение вирусов из организма. При выбранных параметрах оптимального возмущения оно оказалось немонотонной функцией. При этом относительное возмущение было больше в компонентах, отвечающих численности Т-лимфоцитов-эффекторов Ep и предшественников Ee (см. рис. 3 и 4). Начальное возмущение, равное такому оптимальному возмущению, взятому с подходящим коэффициентом ε, вызывает отклик системы, характеризующийся элиминацией вирусов. Величина параметра ε = 0,45 оказывается большой с точки зрения влияния нелинейности на динамику решений, что проявляется в существенном отличии поведения решений задач Коши для линеаризованной системы и исходной нелинейной системы уравнений. Изменение знака параметра ε не оказывает существенного влияния на достижение цели возмущения - элиминацию вирусов до практически нулевых значений, однако уменьшает амплитуду отклика системы (см. рис. 4 и табл. 3). С точки зрения течения вирусной инфекции, такое начальное возмущение состоит из чередующихся обострений и ослаблений инфекционного процесса, которые достигаются за счет комбинации выбранных базисных функций. Для стационарного состояния UII с большой вирусной нагрузкой целью возмущения системы являлось выведение вирусов из организма. На рис. 5 представлены результаты расчетов для сценария, когда в течение 5,6 суток до окончания терапии имело место двукратное введение противовирусных препаратов и имелась возможность 14 раз применить иммуномодулирующую терапию. Оптимальное возмущение оказалось немонотонным по вирусной переменной V и монотонным по переменным Ep и Ee. Взятое с ε = 20 в качестве возмущения начального значения, оно привело к динамике вирусной инфекции, характеризующейся снижением на 4 порядка вирусной нагрузки, активацией иммунного ответа до уровня, удерживающего численность вирусной популяции в области малых значений в течение 200 суток, что сопоставимо с временем жизни экспериментальных животных (мышей). Поведения решений линеаризованной и исходной систем качественно согласуются. Изменение параметров оптимального возмущения привело к качественному изменению его характера - оно стало существенно немонотонным по всем переменным (см. рис. 6). Увеличение в 10 раз веса (вклада) компонент Ep и Ee в максимизируемую норму решения имело следствием уменьшение амплитуды изменения этих компонент решения. При этом, численность вирусной популяции, как и в предыдущем случае (рис. 5), снижалась на 4 порядка, до уровня ниже порога обнаружения стандартными методами, однако существенно быстрее возвращалась к исходному стационарному состоянию. Выравнивание значимости компонент возмущения при одновременном изменении величины и знака параметра ε позволяют, как показано на рис. 7, добиться снижения вирусной нагрузки на 5 порядков, что соответствует практически полному выведению вирусов из организма. При этом, поведение начальных возмущений является существенно немонотонным по всем компонентам, а амплитуда изменений численности клеток иммунной системы остается такой же, как и в предыдущем случае. Из рис. 8 видно, что увеличение числа базисных функций приводит к более немонотонному поведению оптимального возмущения. При этом, снижение вирусной нагрузки происходит в пределах четырех порядков. Увеличение относительного веса компонент Ep и Ee не оказывает влияния на их амплитуды изменения. ТАБЛИЦА 2. Значения параметров модели (3.1), соответствующие исследуемым стационарным состояниям UI и UII. Параметр Размерность UI UII β E0 p bp bd θp θE 1/сутки клетка/мл мл/(частиц·сутки) мл/(частиц·сутки) частицы/мл частицы/мл 1,2 0,08 106 103 7,73 · 10-5 1 · 10-5 7,73 · 10-4 5 · 10-4 3 · 106 10 1 · 105 1,8 · 106 γVE мл/(клетка·сутки) 1,34 · 10-6 Vmvc αEp частицы/мл 1/сутки 4,82 · 107 0,5 αEe 1/сутки 0,1 τ сутки 0,4 τA сутки 5,6 αAP αAE bW (мл/частицы)2/сутки (мл/частицы)2/сутки 1/сутки 7,5 · 10-16 4,36 · 10-14 1 αW 1/сутки 0,11 A B C D РИС. 1. Пример базисных функций для переменных V и W при β = -1 и t0 = -τA (A), β = 1 и t0 = -τA (B), β = -1 и t0 = -τA/2 (C), β = 1 и t0 = -τA/2 (D) A B C D РИС. 2. Пример базисных функций для переменных Ep и Ee при β = -1 и t0 = -τA (A), β = 1 и t0 = -τA (B), β = -1 и t0 = -τA/2 (C), β = 1 и t0 = -τA/2 (D) ТАБЛИЦА 3. Округленные до трех значащих десятичных цифр значения компонент стационарного состояния и максимальные и минимальные значения компонент возмущенного стационарного состояния, соответствующие рис. 3-8. V Ep Ee W Состояние UI 11,5 1,01 · 106 8,96 · 105 104 рис. 3 min 4,85 · 10-36 9,5 · 105 6,5 · 104 6,15 · 10-8 tmin 107,2 0,02 0,04 230 max 1049,4 1,17 · 106 2,39 · 106 2,37 · 103 tmax 5,06 7,63 8,5 6,71 рис. 4 min 3,46 · 10-14 106 3,9 · 105 9,09 · 10-5 tmin 66,42 138,35 151,04 129,36 max 334,74 1,1 · 106 1,73 · 106 1,1 · 103 tmax 158,2 161,67 0,04 160,71 Состояние UII 1,35 · 105 103 5,95 · 104 1,23 · 106 рис. 5 min 23,38 103 108,87 6,93 · 103 tmin 10,33 -4 30,39 56,53 max 5,34 · 105 3,06 · 104 2,38 · 106 1,32 · 106 tmax 0 0 0,79 0,54 рис. 6 min 21,18 656,85 81,7 1,27 · 104 tmin 1,2 · 10-3 10-4 19,75 54,7 max 1,55 · 105 1,06 · 103 6,18 · 104 2,18 · 106 tmax -0,99 -1,56 -1,55 2 · 10-4 рис. 7 min 3,43 975,74 21,51 3,2 · 104 tmin 0,25 10-3 24,32 60,68 max 1,56 · 105 103 6,2 · 104 1,25 · 106 tmax -0,98 -1,5 0 -0,98 рис. 8 min 30,5 693 109 1,56 · 104 tmin 1,2 · 10-3 10-4 19,1 52,9 max 1,56 · 105 1,05 · 103 6,15 · 104 2,18 · 106 tmax -0,81 -1,3 -1,3 2 · 10-4 ТАБЛИЦА 4. Результаты расчета максимальной амплификации. Состояние Номер возмущения topt Γ(topt) Γ�(topt) UI 1 19,0 61,57 2,76 UII 2 5,00 2,82 2,04 UII 3 6,75 44,06 12,82 UII 4 6,75 6,53 4,97 UII 5 7,00 47,78 12,6 11.56 11.54 11.52 11.5 11.48 11.46 5 10.4 10 10.2 10 9.8 9.6 11.44 -6 -4 -2 0 9.4 -6 -4 -2 0 5 12 10 105.5 10 8 105 6 4 104.5 2 0 -6 -4 -2 0 104 -6 -4 -2 0 A 1500 6 1.2 10 1000 1.15 1.1 500 1.05 1 0 0.95 -500 -100 0 100 200 300 0.9 -100 0 100 200 300 6 2.5 10 3000 2 2000 1.5 1 1000 0 0.5 0 -100 0 100 200 300 -1000 -100 0 100 200 300 B РИС. 3. Начальное значение (А) и результаты интегрирования стационарного состояния UI, возмущенного первым оптимальным возмущением с ε = 0,45, темным - интегрирование на основе полных, светлым - на основе линеаризованных уравнений (B) 11.5 6 1.08 10 11.45 1.06 1.04 11.4 1.02 11.35 -6 -4 -2 0 1 -6 -4 -2 0 6 1.8 10 104.5 1.6 1.4 104 1.2 1 103.5 0.8 0.6 -6 -4 -2 0 103 -6 -4 -2 0 A 400 6 1.1 10 300 200 1.05 100 1 0 -100 0.95 -200 -200 0 200 400 600 0.9 -200 0 200 400 600 6 2 10 1500 1.5 1000 500 1 0 0.5 -500 0 -200 0 200 400 600 -1000 -200 0 200 400 600 B РИС. 4. Начальное значение (А) и результаты интегрирования стационарного состояния UI, возмущенного первым оптимальным возмущением с ε = -0,45, темным - интегрирование на основе полных, светлым - на основе линеаризованных уравнений (B) 10 10 8 6 10 5 4 3 2 2 1 0 -6 -4 -2 0 0 -6 -4 -2 0 6 2.5 10 10 6 2 1.5 1 1.227 0.5 0 -6 -4 -2 0 -6 -4 -2 0 A 6 x 10 1 4 x 10 4 0.5 3 0 2 -0.5 1 -1 -100 0 100 200 300 0 -100 0 100 200 300 2.5 6 x 10 6 x 10 2 2 1 1.5 0 1 -1 0.5 -2 0 -3 -0.5 -100 0 100 200 300 -4 -100 0 100 200 300 B РИС. 5. Начальное значение (А) и результаты интегрирования стационарного состояния UII, возмущенного вторым оптимальным возмущением с ε = -20, темным - интегрирование на основе полных, светлым - на основе линеаризованных уравнений (B) 4 6.5 10 6 2.2 10 2 6 1.8 5.5 1.6 1.4 5 1.2 4.5 -6 -4 -2 0 1 -6 -4 -2 0 A 5 x 10 2 1100 1.5 1000 900 1 800 0.5 700 0 -100 0 100 200 300 600 -100 0 100 200 300 4 x 10 8 2.5 6 x 10 6 2 4 1.5 2 1 0 0.5 -2 -100 0 100 200 300 0 -100 0 100 200 300 B РИС. 6. Начальное значение (А) и результаты интегрирования стационарного состояния UII, возмущенного третьим оптимальным возмущением с ε = 0,5371, темным - интегрирование на основе полных, светлым - на основе линеаризованных уравнений (B) 5 2 10 1005 1.5 1000 995 1 990 0.5 985 980 0 -6 -4 -2 0 975 -6 -4 -2 0 4 6.2 10 6.15 6.1 6.05 6 5.95 6 1.3 10 1.25 1.2 1.15 1.1 5.9 -6 -4 -2 0 1.05 -6 -4 -2 0 A 4 x 10 20 15 10 5 0 1005 1000 995 990 985 980 -5 -100 0 100 200 300 975 -100 0 100 200 300 4 x 10 8 5 x 10 15 6 10 4 5 2 0 -100 0 100 200 300 0 -100 0 100 200 300 B РИС. 7. Начальное значение (А) и результаты интегрирования стационарного состояния UII, возмущенного четвертым оптимальным возмущением с ε = -0,3975, темным - интегрирование на основе полных, светлым - на основе линеаризованных уравнений (B) 4 6.5 10 6 2.2 10 2 6 1.8 5.5 1.6 1.4 5 1.2 4.5 -6 -4 -2 0 1 -6 -4 -2 0 A 5 x 10 2 1100 1.5 1000 900 1 800 0.5 700 0 -100 0 100 200 300 600 -100 0 100 200 300 4 x 10 8 2.5 6 x 10 6 2 4 1.5 2 1 0 0.5 -2 -100 0 100 200 300 0 -100 0 100 200 300 B РИС. 8. Начальное значение (А) и результаты интегрирования стационарного состояния UII, возмущенного пятым оптимальным возмущением с ε = 0,4853, темным - интегрирование на основе полных, светлым - на основе линеаризованных уравнений (B) 6. ВЫВОДЫ Современные достижения фармакологии и медицины поставили в повестку исследований разработку подходов к комбинированной терапии неблагоприятных форм течения инфекционных заболеваний. Спектр лечебных воздействий на индивидуальные звенья инфекционного процесса взаимодействия вирусов с организмом хозяина включают противовирусную терапию, иммуномодулирующие воздействия, антифиброзные и гормональные препараты. Поскольку любое инфекционное заболевание представляет собой сложную биологическую систему, динамика которой определяется нелинейными взаимодействиями между различными компонентами инфекции и организма, требуется развитие математических моделей для решения задачи описания и прогнозирования реакции инфекционного процесса на многокомпонентные лечебные воздействия. Построение математических моделей по данным наблюдений является задачей теории идентификации систем. В области математической иммунологии накоплен определенный опыт разработки биологически содержательных моделей инфекционных процессов [8, 33]. Однако возможность построения эффективных мультимодальных воздействий на инфекционный процесс является практически не изученной. Между тем, данная проблема является одной из ключевых для лечения хронических заболеваний в рамках системного подхода, который рассматривает любой патологический процесс как робастную систему по отношению к внешним воздействиям [19, 29]. Например, ВИЧ инфекцию можно рассматривать как робастную систему «вирус-организм хозяина» [16, 28]. Робастность саморегулируемой системы предполагает свойство хрупкости по отношению к определенной комбинации возмущений [7, 19]. Поиск соответствующих возмущений параметров или состояния динамической системы, которые позволили бы разработать эффективные режимы комбинированной терапии, возможен путем построения и анализа свойств чувствительности математической модели заболевания. В данной работе нами предложен алгоритм построения оптимальных возмущений начального состояния математических моделей инфекционных заболеваний, сформулированных в виде систем нелинейных дифференциальных уравнений с запаздывающим аргументом. В качестве базовой модели рассматривалась система из четырех уравнений с двумя запаздываниями, описывающая динамику экспериментальной вирусной инфекции. Комбинированные воздействия моделировались с помощью функций, феноменологически отражающих реакцию системы на применение препаратов с учетом их фармакокинетики. Продемонстрирована принципиальная возможность расчета оптимальных возмущений, приводящих к максимальной амплификации решений по одной или нескольким компонентам инфекционного процесса, например, вирусной нагрузке. Исследована возможность построения оптимальных возмущений для двух типов стационарных состояний, с низким и высоким уровнем вирусной нагрузки [32]. Первый вариант соответствует лечению субпороговых персистентных вирусных инфекций, которые являются проблемой при пересадке органов. Стационарные решения с большим уровнем вирусной нагрузки характерны для инфекций, вызванных вирусами иммунодефицита человека и вирусами гепатита. Поскольку противовирусные и иммуномодулирующие препараты обладают побочным действием, минимизация необходимых доз при сохранении уровня отклика системы является важной частью разработки стратегии лечения. Нами исследована возможность коррекции динамики инфекции и восстановления функции компонент иммунной системы путем малых по норме возмущений стационарных состояний, приводящих к максимальному отклику. Разработанные численные алгоритмы расчета максимального отклика моделей инфекционных заболеваний на мультипараметрические воздействия позволяют перейти к исследованию на регулярной основе принципиально новых подходов направленного воздействия, на основе противовирусных и иммунокорректирующих препаратов, на критические звенья процессов взаимодействия вирусов с организмом человека, и в конечном счете, разработать эффективные методы комбинированной терапии неблагоприятных форм течения инфекционных заболеваний.

Об авторах

Г А Бочаров

Институт вычислительной математики РАН; Российский университет дружбы народов

Email: gbocharov@gmail.com
119333, Москва, ул. Губкина, д. 8; 117198, Москва, ул. Миклухо-Маклая, д. 6

Ю М Нечепуренко

Институт вычислительной математики РАН; Институт прикладной математики им. М. В. Келдыша РАН

Email: yumnech@yandex.ru
125047, Москва, Миусская пл., д. 4

М Ю Христиченко

Институт вычислительной математики РАН; Институт прикладной математики им. М. В. Келдыша РАН

Email: micha.hrist@rambler.ru
125047, Москва, Миусская пл., д. 4

Д С Гребенников

Институт вычислительной математики РАН; Институт прикладной математики им. М. В. Келдыша РАН

Email: dmitry.ew@gmail.com
125047, Москва, Миусская пл., д. 4

Список литературы