Mathematical immunology: processes, models and data assimilation

Full Text

Исследование выполнено за счет гранта Российского научного фонда № 18-11-00171, https://rscf.ru/project/18-11-00171/.

Введение

Иммунная система представляет собой сложный многомасштабный и мультифизический объект, изучение закономерностей функционирования которого, в рамках системного подхода, предполагает активное использование методов математического моделирования. Разработка математических моделей иммунологических процессов связана с решением следующих задач: (1) построение системы уравнений, (2) выбор содержательного критерия близости модели и данных наблюдений, (3) идентификация параметров и оценка неопределенности параметров, (4) выбор оптимальной модели и (5) исследование локальной и глобальной чувствительности. В данной работе излагаются перспективные подходы в этой области, связанные с моделированием нового класса мультифизических процессов иммунной системы: миграция клеток в лимфатических узлах, лимфодинамика, гомеостатическая регуляция иммунного ответа при хронических инфекциях. В силу различия природы данных процессов, для их описания применяются различные классы математических моделей, при построении которых используются системы обыкновенных дифференциальных уравнений (ОДУ), стохастических дифференциальных уравнений (СДУ), клеточные автоматы (КА) и нейронные сети (НС). Активное использование данных технологий моделирования позволит существенно расширить спектр вычислительных инструментов для решения задач моделирования и анализа системной иммунологии.

Материалы и методы

Разработка математической модели основана на выделении механизменно- значимых компонент иммунных процессов (переменные пространства состояний системы), формирования причинно-следственных взаимоотношений между ними в виде биологической схемы взаимодействий и переходов между состояниями [3]. Преобразование схемы в систему уравнений математической модели характеризуется множественностью возможных реализаций, и желательно сузить класс возможных решений. Наличие экспериментальных и клинических данных существенно влияет на объективизацию этого перехода, позволяя найти оптимальный баланс между сложностью модели и информационным содержанием блока эмпирических данных, которые необходимо описать или усвоить с помощью модели. Усвоение данных наблюдений в математической модели относится к процедуре идентификации или решения обратной задачи. Задачи идентификации делятся на коэффициентные обратные задачи, когда параметры модели определяются по значениям переменных, и проблемы структурной идентификации, каждая из которых реализуется в рамках определенной методологии.

Метод максимального правдоподобия (ММП)

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

Байесовский метод (БП)

В отличие от ММП, в рамках которого неизвестный параметр считается детерминированной величиной, в Байесовском подходе искомый параметр модели считается случайной величиной с некоторым законом распределения. При этом для нахождения Байесовской оценки используется дополнительная априорная информация. Сочетание априорной информации о параметрах и информации о данных наблюдений позволяет получить апостериорную оценку их закона распределения.

Матрица Фишера

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

Критерии сложности моделей

Использование ММП или БП предполагает фиксированную структуру модели, т. е. способ параметризации в виде уравнений описываемых моделью процессов. В случае наличия разных гипотез о свойствах процессов иммунной системы можно сформировать множество соответствующих им моделей. Для выбора оптимальной по степени сложности модели существуют информационные критерии, такие как критерий Акаике и критерий минимальной длины описания. Оптимальной моделью является наиболее простая по сложности модель, которая тем не менее обладает хорошей степенью согласия с экспериментальными данными. Для использования критерия минимальной длины описания требуется вычисление матрицы Фишера [6].

Результаты и обсуждение

Миграция клеток – многокритериальный подход

Для описания пространственно-временной динамики развития иммунного ответа в лимфатических узлах (ЛУ), требуется построение математических моделей подвижности и миграции клеток. В работе [4] была реализована и откалибрована математическая модель движения лимфоцитов на основе II закона Ньютона. На каждую клетку действуют три вида сил: (1) стохастическая сила активной подвижности клетки, неявно описывающая взаимодействия с ретикулярными структурами лимфоидной ткани, (2) силы специфических и неспецифических взаимодействий с другими клетки, параметризованная сила адгезии между мембранами клеток, (3) диссипативная сила трения в вязкой среде, пропорциональная значению мгновенной скорости клетки. Для калибровки модели необходимо выбрать критерий согласия с экспериментальными данными. Для описания траекторий клеток были использованы доступные в литературе эмпирические распределения трех характеристик движения лимфоцитов в ЛУ: поступательных скоростей клеток, скоростей углов поворота клеток, индексов меандрирования клеток (рис. 1). Для набора статистики значения характеристик рассчитывались раз в 30 секунд в течение продолжительного времени (5 часов). Для оценки степени отклонения между экспериментальными и предсказанными моделью характеристиками использовалось расстояние Колмогорова–Смирнова между соответствующими эмпирическими функциями распределениями.

 

Рисунок 1. Калибровка модели движения иммунных клеток в лимфоузле по трем характеристикам траекторий: поступательным скоростям, скоростям углов поворота и индексов меандрирования

Figure 1. Calibration of the model of immune cell motility in lymph node by the three characteristics of cell trajectories: translation speeds, turning angle speeds, meandering indices

 

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

Лимфоток через ЛУ – обучение нейронной сети

Расчет задачи течения лимфы в ЛУ является сложным с вычислительной точки зрения процессом. Учитывая разнообразие размеров, форм и структур внутренней организации, а также вариативность граничных условий, таких как поступающий в лимфатический узел поток лимфы и баланс давлений на выходе из него и во внутренней кровеносной системе, решение задач течения лимфы в ЛУ становится существенным препятствием для проведения расчетов лимфодинамики и связанных с ней процессов в организме человека. В работе [7] предложен подход к моделированию лимфотока в ЛУ с замещением полноценной расчетной физической модели, требующей построения сеточной аппроксимации лимфатического узла и решения системы уравнений, на искусственную нейронную сеть (ИНС), обученную на наборе заранее сформированных результатов вычислений исходной модели. Зафиксировав форму лимфатического узла, определяющими результаты расчетов можно считать всего лишь 6 параметров: L – скорость абсорбции лимфы в капиллярах; α₁ – гидравлическое сопротивление в субкапсулярном синусе; α₂ – гидравлическое сопротивление в кортексе и медуллярной зоне лимфатического узла; Q₁ – афферентный (входящий) поток лимфы; P₂ – давление в эфферентном (выходящем) сосуде; P_ν – давление в кровеносных капиллярах (рис. 2).

 

Рисунок 2. Схема модели на основе искусственной нейронной сети и процедура идентификации ее параметров путем минимизации отклонений предсказания нейросети X_NN от решения X_IE механизменной модели

Примечание. На рисунке слева приведены: схема нейросети, система реализующих ее уравнений (x_h – выходные значения скрытого слоя, W₁, W₂ – матрицы весовых коэффициентов для скрытого и результирующего слоя, b₁, b₂ – свободные коэффициенты для скрытого и результирующего слоя, a₁, a₂ – нормализующие коэффициенты для входных значений, ξ(x) – нелинейная функция активации); справа представлены принцип обучения нейронной сети методом обратного распространения ошибки, формула оценки ошибки результата вычислений.

Figure 2. ANN-based model scheme and it's parameter identification procedure by minimization of the difference between network prediction X_NN from initial mechanistic model solution X_IE

Note. Figure on the left shows: a topological scheme of a neural network, a system of equations that implements it (x_h, output values of the hidden layer; W₁, W₂, matrices of weight coefficients for the hidden and resulting layers; b₁, b₂, free coefficients for the hidden and resulting layers; a₁, a₂, normalizing coefficients for input values; ξ(x), non-linear activation function); on the right are the principle of training the neural network by the method of error backpropagation, the formula for estimating the error of the calculation result.

 

Используя их как входные параметры, была сформирована схема искусственной нейронной сети, содержащей: 6 нейронов входного слоя (по числу определяющих результат параметров), 6 нейронов скрытого слоя с нелинейной функцией активации (сигмоидальная неотрицательная) и двух нейронов выходного слоя, с которых снимались результаты вычислений физической модели – поток в эфферентном сосуде и давление в афферентном. При этом использовался классический принцип построения связей между нейронами: каждый нейрон последующего слоя имел связи со всеми нейронами слоя предыдущего плюс свободный коэффициент. Величина функции ошибки, т. е. невязка между решением физической и нейросетевой модели, оценивалась как относительное отклонение полученного решения от целевого.

Результаты, полученные при тестировании обученной сети на данных, не принадлежащих обучающей выборке, показали, что при несущественном ухудшении качества результатов, использование нейронной сети позволяет многократно ускорить последующие расчеты: вместо 4-4,5 минуты на решение методом интегральных уравнений, ИНС позволяет получить результат всего за несколько микросекунд, что дает возможность использовать подобные методы приближения результатов в составных проектах, включающих дренаж лимфатических узлов как один из элементов более комплексных моделей (ЛС человека).

Противовирусный иммунный ответ – приближение данных кинетики и допустимые целевые множества

Для решения обратной задачи идентификации параметров модели по экспериментальным данным часто требуется задание, в том или ином виде, дополнительных ограничений. В работе [5] была проведена калибровка модели противовирусного иммунного ответа Марчука–Петрова для описания динамики новой коронавирусной инфекции COVID-19. Модель представляет собой систему ОДУ с запаздывающими аргументами из 13 переменных, описывающую вирусную нагрузку, врожденный интерфероновый ответ, развитие Т-клеточного и гуморального специфического иммунного ответа. В модели параметризовано усиление специфического иммунного ответа за счет воспаления и его подавление за счет разрушения клеток легких. В качестве экспериментальных данных для калибровки использовались опубликованные данные кинетики вирусной нагрузки в носоглотке инфицированных добровольцев, а также дополнительные литературные данные о диапазоне значений концентрации интерферона, антител и цитотоксических лимфоцитов в крови в некоторые моменты времени. Для оценки параметров модели решалась задача минимизации функционала относительных отклонений вирусной нагрузки Ф. Для всех параметров модели были оценены диапазоны допустимых значений параметров на основании литературных данных о скоростях иммунных процессов. Значения большинства параметров модели были зафиксированы равными оценкам параметров модели, откалиброванной ранее для описания динамики инфекции вирусами гриппа А, либо равными оценкам на основе литературных данных. Для оставшихся параметров был проведен локальный анализ чувствительности функционала Ф и его частей, учитывающих отклонения вирусной нагрузки на разных временных интервалах. Были выделены параметры, наиболее влияющие на разные стадии развития инфекции. Для калибровки решалась серия задач оптимизации с двухсторонними ограничениями для подгрупп параметров, влияющих на Ф последовательно на каждом временном интервале. Для валидации корректности получаемых оценок параметров проводилось сравнение решений модели с экспериментальными данным по другим переменным модели.

Т-клеточный гомеостаз – регуляции баланса деления и дифференцировки

Отличительным признаком хронических вирусных инфекций и онкологических заболеваний являются нарушения работы иммунных механизмов, в частности, истощение T-лимфоцитов. Истощенные Т-лимфоциты экспрессируют ингибиторные рецепторы (такие как PD1, TIM3), обладают сниженной способностью к пролиферации, выделяют меньше цитокинов. В недавней работе [2] были выделены 4 субпопуляции истощенных Т-лимфоцитов, которые соответствуют различным стадиям дифференцировки и обладают различными характеристиками: T_PROG1 (LY108⁺CD69⁺) – покоящиеся, резидентные, T_PROG2: (LY108⁺ CD69^-) – циркулирующие, пролиферирующие, T_EXINT: (LY108^-CD69^-) – циркулирующие, обладающие слабым цитотоксическим эффектом, T_EXTER: (LY108^-CD69⁺) – покоящиеся, терминально истощенные. Баланс деления и дифференцировки данных субпопуляций поддерживается, в частности, за счет взаимодействия с дендритными клетками (DC). Так, SIRPa⁺ PD-L1⁺ DC ингибируют пролиферацию, а XCR⁺1 DC частично восстанавливают функциональность истощенных Т-лимфоцитов, способствуют делению их вместо дифференцировки. Истощение Т-лимфоцитов является частично обратимым за счет применения антител, блокирующих взаимодействие ингибиторных рецепторов и их лигандов. Предполагаемая биологическая схема регуляции дифференцировки и деления может быть сформулирована в виде системы ОДУ (рис. 3). Построенную математическую модель можно использовать, например, для анализа экспериментальных данных, описывающих гомеостаз системы. Для этого требуется идентифицировать параметры системы ОДУ, соответствующие ее стационарному состоянию. Описанный подход можно применить для оценки влияния терапии (например, aPD-L1) на параметры модели, характеризующие те или иные биологические процессы. Применение математического моделирования для прогнозирования эффекта иммунотерапии при ВИЧ-инфекции описано, например, в [8]. В случае анализа стационарных состояний системы нам требуется найти значения параметров, при которых правая часть системы равна нулю. В случае, если построенная модель является линейной, для решения задачи может применяться метод линейного программирования. В отличие от решения системы линейных алгебраических уравнений, при использовании этого подхода решается задача оптимизации с условиями, что, с одной стороны, позволяет учесть ограничения на параметры, с другой – искать вектор с минимальной нормой. Схема математического моделирования эффекта применения иммунотерапии на основе идентификации параметров Т-клеточного гомеостаза приведена на рисунке 3.

 

Рисунок 3. Применение математического моделирования для оценки влияния иммунотерапии на параметры Т-клеточного гомеостаза

Figure 3. Application of mathematical modelling for estimation of immunotherapy effect on T cell homeostasis parameters

 

Заключение

Представлены современные подходы к разработке математических моделей, количественно согласованных с эмпирическими данными и характеристиками иммунных процессов различной природы.

About the authors

Dmitry S. Grebennikov

G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; I. Sechenov First Moscow State Medical University

Email: dmitry.ew@gmail.com

Junior Research Associate, G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; Research Associate, I. Sechenov First Moscow State Medical University, Moscow, Russian Federation

Russian Federation, Moscow; Moscow

Valeriya V. Zheltkova

G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; I. Sechenov First Moscow State Medical University

Author for correspondence.
Email: zheltkova_v_v@staff.sechenov.ru

Junior Research Associate, G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; Junior Research Associate, I. Sechenov First Moscow State Medical University, Moscow, Russian Federation

Russian Federation, Moscow; Moscow

Rostislav S. Savinkov

G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; I. Sechenov First Moscow State Medical University

Email: savinkov_r_s@staff.sechenov.ru

Junior Research Associate, G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; Junior Research Associate, I. Sechenov First Moscow State Medical University, Moscow, Russian Federation

Russian Federation, Moscow; Moscow

Gennady A. Bocharov

G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; I. Sechenov First Moscow State Medical University

Email: gbocharov@gmail.com

PhD, MD (Phys.-Math.), Leading Research Associate, G. Marchuk Institute of Numerical Mathematics, Russian Academy of Sciences; Professor, I. Sechenov First Moscow State Medical University, Moscow, Russian Federation

Russian Federation, Moscow; Moscow

References

Supplementary files