1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 60, № 11, 2020

Журнал вычислительной математики и математической физики, 2020, T. 60, № 11, стр. 1950-1961

Математическое моделирование эпидемии Уханьского коронавируса COVID-2019 и обратные задачи

С. И. Кабанихин 1, 2, О. И. Криворотько 1, 2, *

1 Институт вычислительной математики и математической геофизики СО РАН
630090 Новосибирск, пр-т Акад. Лаврентьева, 6, Россия

2 Новосибирский государственный университет
630090 Новосибирск, ул. Пирогова, 2, Россия

* E-mail: krivorotko.olya@mail.ru

Поступила в редакцию 02.03.2020
После доработки 02.03.2020
Принята к публикации 07.07.2020

Полный текст (PDF)

Аннотация

В работе рассмотрены математические модели распространения нового Уханьского коронавируса COVID-2019, вспышка эпидемии которого началась в декабре 2019 г. в провинции Ухань. Для контроля эпидемиологической ситуации необходимо разработать соответствующие математические модели. В работе приведен обзор математических моделей распространения COVID-2019, описываемых системами нелинейных обыкновенных дифференциальных уравнений (ОДУ). Некоторые коэффициенты и начальные данные систем ОДУ неизвестны или заданы усредненно. Задача идентификации параметров моделей сводится к минимизации квадратичного целевого функционала. В силу нелинейности ОДУ решение обратных задач эпидемиологии может быть неединственно, поэтому в работе изложены подходы к исследованию идентифицируемости обратных задач, которые позволяют установить, какие из неизвестных параметров (либо их комбинаций) могут быть однозначно и устойчиво восстановлены по имеющейся дополнительной информации. Приведены методы решения задачи минимизации, основанные на комбинации глобальных (методы покрытий, природоподобные алгоритмы, многоуровневые градиентные методы) и локальных методов (градиентные методы, метод Нелдера-Мида). Библ. 44. Фиг. 7. Табл. 2.

Ключевые слова: математические модели, COVID-2019, коронавирус, эпидемиология, обратные задачи, оптимизация, регуляризация, идентифицируемость, ОДУ, тензорное разложение, природоподобные алгоритмы, градиентные методы.

1. ВВЕДЕНИЕ

В декабре 2019 г. произошла вспышка пневмонии в Ухане 2019–2020 гг., в результате которой был впервые обнаружен штамм COVID-2019 при анализе нуклеиновой кислоты у пациента с пневмонией. Репродуктивное число по данным Китайского центра по контролю и профилактике заболеваний оценивается между 2 и 3, что соответствует количеству людей, которые заражаются от одного инфицированного, и пока оно больше 1, эпидемия будет расти (см. фиг. 1, 2).

Фиг. 1.

Количество инфицированных COVID-2019 (синяя линия), вылеченных (оранжевая линия) и умерших (зеленая линия) индивидуумов с 22.01.2020 по 11.02.2020 на территории Китая.

Фиг. 2.

Зарегистрированные случаи обнаружения Уханьского коронавируса COVID-2019 24 января 2020 г. (слева) и 10 февраля 2020 г. (справа) на территории Китая [3], [7].

Возможным источником вируса COVID-2019 являются летучие мыши, поскольку РНК образцов COVID-2019 на 96% совпала с РНК вируса, который ранее находили у азиатских подковоносов [1]. Репликация вируса происходит преимущественно в нижних дыхательных путях, вызывая большой выброс цитокинов и иммунный ответ в организме, при котором наблюдается снижение количества Т-лимфоцитов в крови, отвечающих за защитные функции организма [2]. По состоянию на 11 февраля 2020 г. число подтвержденных инфицированных случаев составило 43 143 человека, из которых зарегистрировано 1018 летальных исходов и 4347 вылеченных людей [3]. Генетически вирус COVID-2019 на 80% идентичен тяжелому острому респираторному синдрому (SARS), вспышка которого наблюдалась в Китае в 2003 г., в результате которой было инфицировано более 5000 человек [4] (см. табл. 1). Однако скорость распространения коронавируса COVID-2019 значительно выше, чем в случае SARS, поскольку с 24 января количество зарегистрированных случаев в провинции Ухань выросло за 17 дней с 549 человек до 31  728 человек (см. фиг. 2).

Таблица 1.  

Статистическая информация о распространенности коронавирусов SARS (на 31 июля 2003 г.) [4] и COVID-2019 (на 10 февраля 2020 г.) [7]

Вирус Инфицированные (чел.) Вылеченные (чел.) Смертность (чел.) Инкубационный период (дни)
SARS 5327 – 349 0
COVID-2019 43 143 4347 1018 5.2

Также ключевым эпидемиологическим параметром является инкубационный период, который оценивается в 3–7 дней по оценкам ВОЗ. В работах [5], [6] получено оптимальное значение инкубационного периода 5.2 дня.

30 января Всемирная организация здравоохранения объявила вспышку эпидемии, связанную с COVID-2019, чрезвычайной ситуацией в области здравоохранения международного значения [7].

В связи со сложившейся ситуацией необходимо прогнозировать развитие эпидемиологической ситуации в Китае и в мире. Одним из способов предсказания распространения коронавируса COVID-2019 является математическое моделирование распространения инфекции в популяции с заданным источником инфекции, скоростями перехода между группами людей со схожими характеристиками (симптомы, карантин, антитела, штаммы), смертности, латентного периода, степени изоляции (миграция населения между провинциями, запрет на ввоз продуктов) и статистической информации о зараженных и вылеченных индивидуумах в фиксированные моменты времени. Математическая модель, основанная на законе баланса масс и описываемая системой нелинейных обыкновенных дифференциальных уравнений (ОДУ), наиболее точно описывает распространение инфекционных заболеваний в популяции, разделенной на группы индивидуумов со схожими характеристиками (например, чувствительные, инфицированные, находящиеся на лечении, вылеченные и т.п.). Большинство коэффициентов и начальных условий ОДУ не известно или может быть оценено с большой погрешностью. Например, в работе [8] приведен прогноз количества зарегистрированных случаев вируса COVID-2019 на основании данных на период 23 января 2020 г., который составил 75 000 человек на 10 февраля 2020 г., хотя реальные данные – 40 645 человек. Такая погрешность прогнозирования развития инфекционного заболевания влечет к огромным затратам принятия мер по контролю эпидемиологической ситуации, а также к ошибочным выводам. Для идентификации неизвестных параметров математических моделей систем ОДУ необходимо использовать дополнительную информацию о количестве зараженных и вылеченных индивидуумов по каждому дню. Используемый в работе подход является классическим в обратных задачах [9]. В следующих разделах мы опишем существующие математические модели распространения коронавируса COVID-2019, приведем методы решения возникающих обратных задач и примеры решения.

Статья организована следующим образом. В разд. 2 приведены математические модели эпидемиологии COVID-2019 и постановки обратных задач. Разд. 3 посвящен описанию идентифицируемости математической модели. В разд. 4 приведена классификация методов оптимизации решения обратных задач и область их применения.

2. МАТЕМАТИЧЕСКИЕ МОДЕЛИ ДИНАМИКИ COVID-2019

Математические модели эпидемиологии основаны на камерной структуре, описываются системами нелинейных ОДУ

(1)
$\dot {X} = F\left( {X(t),\varphi } \right),~~~X({{t}_{0}}) = {{X}_{0}},\quad t \in \left( {{{t}_{0}},T} \right),$
и характеризуют изолированное распространение инфекционного заболевания внутри региона. Вспышка очередного инфекционного заболевания влечет разработку новых математических моделей эпидемиологии, учитывая особенности течения заболевания. Одна из первых моделей эпидемии туберкулеза, построенная в работе [10] на основе опыта использования математических моделей в метеорологии, демографии, экономике и эпидемиологии острых заболеваний, “схватывает” главное отличие эпидемического процесса туберкулеза от эпидемических процессов острых болезней – наличие длительного латентного периода. Далее развитие модели Ваалера было получено в работах советских математиков Г.И. Марчука и его последователей [11], [12]. В частности, широкое применение получила математическая модель распространения пневмонии [13].

В работе [1] разработана математическая модель распространения коронавируса COVID-2019 от предполагаемого источника инфекции (летучих мышей) к человеку. Поскольку отследить источник инфекции до сих пор не удалось, авторы адаптировали математическую модель передачи вируса от рынка морепродуктов людям. Модель, схематично изображенная на схеме (см. фиг. 3), состоит из 14 ОДУ типа (1) и характеризуется 25 коэффициентами φ, описывающими скорости перехода между 14 группами, разделение которых основано на известной SEIR модели, описывающей динамику инфекции между восприимчивыми, подверженными заражению, инфицированными и вылеченными индивидуумами. Задача идентификации неизвестных коэффициентов по доступным статистическим данным о количестве инфицированных и вылеченных индивидуумов в фиксированные дни некорректна, то есть ее решение не единственно и не устойчиво [14], [15]. При использовании такой модели с целью прогнозирования необходимо получить больше статистических данных.

Фиг. 3.

Блок-схема математической модели источник–рынок–люди [1].

Более подробная математическая модель метапопуляции основана на глобальной сети локальных групп населения, соединенных ребрами, представляющими пассажиропотоки между городами [16]. В каждом узле сети авторы локально моделируют динамику вспышек с использованием камерной SEIR модели с дискретным временем, близкую по своей структуре к математической модели эпидемии туберкулеза 1962 г. [10]. Параметры SEIR оцениваются на основе 5-дневного периода инкубации. Предполагается, что первоначальные случаи COVID-2019 присутствуют только в Ухани, и пограничный контроль не учитывается.

В работе [17] авторы предложили детерминистическую SEIR камерную модель динамики нового коронавируса с целью оценить влияние вмешательств общественного здравоохранения на инфекцию (см. фиг. 4). Согласно модели, вспышка должна была достичь максимума 7 февраля 2020 г. со значительным низким пиковым значением. Предложенная эпидемиологическая модель SEIR основана на 8 ОДУ с 12 коэффициентами φ, и задача подбора параметров решена с помощью метода интервальной генерации [18]. Отметим, что прогноз был получен с погрешностью, так как пиковое значение в рамках роста еще не получено на 12 февраля 2020 г., что требует разработки более точных методов идентификации параметров систем ОДУ.

Фиг. 4.

Диаграмма математической модели распространения COVID-2019 [17].

Новая математическая модель интегро-дифференциальных уравнений с запаздыванием, описывающая эпидемию коронавируса COVID-2019 в Китае была предложена в работе [19]. Модель состоит из четырех интегро-дифференциальных уравнений, в которых определяются три неизвестных коэффициента локальным градиентным методом по дополнительной информации о количестве инфицированных и вылеченных индивидуумов за две недели. В силу неединственности решения задачи идентификации, полученные результаты не гарантируют достоверность расчетов ввиду применения локальных итерационных методов [20], [21].

Тот же подход построения камерных математических моделей, основанных на системах нелинейных ОДУ, был предложен для других близких острых инфекционных заболеваний, таких как вспышка SARS в Китае 2003 г. [22]–[24], ближневосточного респираторного синдрома (MERS) в Саудовской Аравии и Южной Корее [25], туберкулеза и ВИЧ в России и Казахстане [26]–[28] и др. Несмотря на геномную схожесть COVID-2019 с коронавирусом SARS и MERS, течение заболевания, включая инкубационный период, сильно разнится, что делает невозможным применение разработанных моделей к новому коронавирусу COVID-2019. Наличие латентного периода у COVID-2019 схоже с математическими моделями динамики туберкулеза, на основании которых авторы настоящей статьи получили верифицированный на статистических данных прогноз [26]–[28].

2.1. Постановка обратной задачи

Учитывая меры, принятые в Китае по сдерживанию эпидемии, включая карантин и изоляцию, разделим популяцию в выбранной провинции на восприимчивых (S), подверженных заражению (E), инфицированных предсимптомных (не симптоматические) (A) и с симптомами (I), госпитализированных (H) и вылеченных индивидуумов (R). Также выделим отдельно группы на карантине: восприимчивые (Sq) и подверженные заражению (Eq) [17]. Диаграмма камерной математической модели представлена на фиг. 4. Учитывая закон баланса масс, математическая модель для системы ОДУ (1) запишется в виде

$\begin{gathered} \dot {S} = - \left( {\beta c + cq\left( {1 - \beta } \right)} \right)S\left( {I + \theta A} \right) + \lambda {{S}_{q}}, \\ \dot {E} = \beta c\left( {1 - q} \right)S\left( {I + \theta A} \right) - \sigma E, \\ \dot {I} = \sigma \rho E - \left( {{{\delta }_{I}} + \alpha + {{\gamma }_{I}}} \right)I, \\ \end{gathered} $
(2)
$\begin{gathered} \dot {A} = \sigma \left( {1 - \rho } \right)E - {{\gamma }_{A}}A, \\ {{{\dot {S}}}_{q}} = \left( {1 - \beta } \right)cqS\left( {I + \theta A} \right) - \lambda {{S}_{q}}, \\ \end{gathered} $
$\begin{gathered} {{{\dot {E}}}_{q}} = \beta cqS\left( {I + \theta A} \right) - {{\delta }_{q}}{{E}_{q}}, \\ \dot {H} = {{\delta }_{I}}I + {{\delta }_{q}}{{E}_{q}} - \left( {\alpha + {{\gamma }_{H}}} \right)H, \\ \dot {R} = {{\gamma }_{I}}I + {{\gamma }_{A}}A + {{\gamma }_{H}}H, \\ \end{gathered} $
с начальными данными

(3)
$\begin{gathered} S({{t}_{0}}) = {{S}_{0}},\quad E({{t}_{0}}) = {{E}_{0}},\quad I({{t}_{0}}) = {{I}_{0}},\quad A({{t}_{0}}) = {{A}_{0}}, \\ {{S}_{q}}({{t}_{0}}) = {{S}_{{{{q}_{0}}}}},\quad {{E}_{q}}({{t}_{0}}) = {{E}_{{{{q}_{0}}}}},\quad H({{t}_{0}}) = {{H}_{0}},\quad R({{t}_{0}}) = {{R}_{0}}. \\ \end{gathered} $

Описание параметров математической модели (2) и значения данных (3) в случае провинции Ухань на 10 января 2020 г. приведены в табл. 2.

Таблица 2.  

Описание параметров и их приближенных значений, а также начальных данных для математической модели распространения коронавируса COVID-2019 в Китае (2), (3) [17]

Параметр Описание Среднее значение
$c$ Частота контактов 14.781
$\beta $ Вероятность передачи инфекции индивидууму $2.1011 \cdot {{10}^{{ - 8}}}$
$q$ Карантинная частота заражения индивидуумов $1.8887 \cdot {{10}^{{ - 7}}}$
$\sigma $ Коэффициент перехода зараженных лиц к инфицированному классу 1/7
$\lambda $ Частота, с которой помещенные на карантин незараженные контакты были переданы в сообщество 1/14
$\rho $ Вероятность наличия симптомов среди инфицированных лиц 0.86834
${{\delta }_{I}}$ Коэффициент перехода зараженных симптомами индивидов к инфицированному карантинному классу 0.13266
${{\delta }_{q}}$ Скорость перехода подвергнутых карантину лиц к карантинному зараженному классу 0.1259
${{\gamma }_{I}}$ Скорость выздоровления симптомно инфицированных лиц 0.33029
${{\gamma }_{A}}$ Скорость выздоровления бессимптомно инфицированных лиц 0.13978
${{\gamma }_{H}}$ Коэффициент выздоровления инфицированных на карантине лиц 0.11624
$\alpha $ Уровень смертности от заболевания $1.7826 \cdot {{10}^{{ - 5}}}$
Начальные значения Описание Средние значения
${{S}_{0}}$ Начальная восприимчивая популяция провинции Ухань 11 081 000
${{E}_{0}}$ Первоначально выявленное население 105.1
${{I}_{0}}$ Начальная симптоматическая инфицированная популяция 27.679
${{A}_{0}}$ Начальная бессимптомная инфицированная популяция 53.839
${{S}_{{{{q}_{0}}}}}$ Первоначальная восприимчивая популяция на карантине 739
${{E}_{{{{q}_{0}}}}}$ Первоначально помещенное в карантин подверженное заражению население 1.1642
${{H}_{0}}$ Первоначально инфицированное население на карантине 1
${{R}_{0}}$ Первоначально вылеченные индивидуумы 2

В математической модели (2), (3) присутствует 14 неизвестных параметров $p = (c,\beta ,q,\rho ,{{\delta }_{I}},{{\delta }_{q}},$ ${{\gamma }_{I}},{{\gamma }_{A}},{{\gamma }_{H}},\alpha ,~{{E}_{0}},{{I}_{0}},{{A}_{0}},{{E}_{{{{q}_{0}}}}}) \in {{\mathbb{R}}^{{14}}}.$ Для идентификации указанных параметров сформулируем обратную задачу, которая будет состоять в определении вектора параметров p по доступным дополнительным измерениям о количестве восприимчивых в фиксированные дни tk:

(4)
$S({{t}_{k}}) = {{S}_{k}},\quad {{S}_{q}}({{t}_{k}}) = {{S}_{{{{q}_{k}}}}},\quad H({{t}_{k}}) = {{H}_{k}},\quad R({{t}_{k}}) = {{R}_{k}},\quad {{t}_{k}} \in \left( {{{t}_{0}},T} \right),\quad k = 1, \ldots ,K.$

В случае векторной записи (1): $X = (S,E,I,A,{{S}_{q}},{{E}_{q}},H,R)$, $p = (\varphi ,~{{E}_{0}},{{I}_{0}},{{A}_{0}},{{E}_{{{{q}_{0}}}}})$. Запишем данные (4) в следующем виде:

${{X}_{i}}({{t}_{k}}) = {{f}_{{ik}}},\quad k = 1, \ldots ,K,\quad i \subset \{ 1, \ldots ,N\} ,$
где ${{f}_{{ik}}} = ({{S}_{k}},{{S}_{{{{q}_{k}}}}},{{H}_{k}},{{R}_{k}})$, $N = 8$.

В операторной форме обратная задача (2)–(4) запишется в виде

$A(p) = f,$
$A{\kern 1pt} :Q \to F$, $Q = \{ p \geqslant 0\} $, F – гильбертовы пространства допустимых параметров и измерений, соответственно, $f = ({{f}_{{i1}}},{{f}_{{i2}}}, \ldots ,{{f}_{{iK}}})$ – вектор данных, $i = 1,5,7,8$.

Обратная задача (2)–(4) в общем случае некорректна, то есть ее решение может быть не единственным или неустойчивым (производная Фреше оператора обратной задачи $A{\kern 1pt} '$ не имеет обратного). В разд. 3 приведен анализ идентифицируемости, который позволяет выявить набор идентифицируемых параметров по имеющимся данным (4) для дальнейшего их определения [15], [14], [29].

Обратная задача (2)–(4) сводится к задаче минимизации целевого функционала $J(p) = \left\langle {A(p) - f,A(p) - f} \right\rangle $, который в нашем случае имеет вид

$J(p) = \frac{1}{K}\mathop \sum \limits_{k = 1}^K \,\mathop \sum \limits_i \,{{({{X}_{i}}({{t}_{k}};p) - {{f}_{{ik}}})}^{2}},$
и характеризует квадратичное отклонение модельных данных от статистических. Таким образом, обратная задача (2)–(4) записывается в следующем виде:

(5)
$p{\kern 1pt} * = \mathop {\operatorname{argmin} }\limits_{p \in Q} J(p).$

В следующем разделе приводится классификация методов решения задачи (5).

3. ИДЕНТИФИЦИРУЕМОСТЬ МАТЕМАТИЧЕСКОЙ МОДЕЛИ

В некоторых случаях решение обратной задачи (2)–(4) может быть неединственно или неустойчиво по отношению к ошибкам в измерениях (4). С целью определения единственного набора параметров (или их комбинаций), а также степени чувствительности неизвестных параметров к ошибкам в данных обратной задачи проводится анализ идентифицируемости математических моделей для систем ОДУ [14], [30–31]. Условно разделяют три вида идентифицируемости (диаграмму см. на фиг. 5):

Фиг. 5.

Классификация методов идентифицируемости [14].

1. Априорная (структурная) идентифицируемость используется при изучении структуры модели в случае непрерывных данных. К методам структурной идентифицируемости относятся методы придаточной функции, дифференциальной алгебры, разложения в ряд, теории графов и другие [14], [30]. Например, метод, основанный на теории графов, не только определяет идентифицируемость модели, но и позволяет найти замену переменных специального вида, приводящую исходную модель к идентифицируемой.

2. Практическая (апостериорная) идентифицируемость позволяет выявить идентифицируемые параметры относительно зашумленных данных методами Монте-Карло, корреляционной матрицы и другими [14], [29].

3. Анализ чувствительности выявляет последовательность чувствительных параметров к ошибкам в измерениях и позволяет оценить степень чувствительности [14], [29].

Для исследования идентифицируемости обратной задачи (2)–(4) мы используем программное обеспечение DAISY [31]. В результате анализа задача идентификации 14 параметров $p$ по измерениям (4) имеет неединственное решение, то есть модель локально неидентифицируема. При сокращении вектора параметров $(\beta ,q,\rho ,{{\delta }_{I}},{{\delta }_{q}},{{E}_{0}},{{I}_{0}},{{A}_{0}},{{E}_{{{{q}_{0}}}}}) \in {{\mathbb{R}}^{9}}$ математическая модель (2)–(4) становится локально идентифицируемой.

4. МЕТОДЫ ОПТИМИЗАЦИИ

Методы оптимизации условно можно разделить на глобальные, локальные и комбинированные. К первым относятся природоподобные алгоритмы, включающие генетический алгоритм, метод имитации отжига, метод роя частиц, дифференциальной эволюции и т.п., а также методы покрытий, такие как метод неравномерных покрытий Ю.Г. Евтушенко [32]–[34], тензорное разложение. Глобальные методы оптимизации позволяют локализовать область экстремума и эффективны в пространствах большой размерности. Локальные методы оптимизации подразделяются на две группы: нулевого порядка, использующие только значение самого функционала (методы Нелдера–Мида, Хука–Джейвса) и градиентные методы (Гаусса–Ньютона, Левенберга–Маквардта и т.п.). Особенность локальных методов состоит в высокой производительности и сходимости к локальному экстремуму для градиентных методов. Главным недостатком локальных методов является сходимость к локальному минимуму, особенно в пространствах большой размерности. Комбинированные методы основаны на следующей идее: методы глобальной оптимизации определяют область глобального экстремума, стартуя в которой локальные методы уточняют сам экстремум [26]. Диаграмма условной классификации методов оптимизации представлена на фиг. 6. Далее приведем основные методы трех групп: глобальные стохастические (п. 3.1), методы покрытий (п. 3.2) и локальные методы градиентного типа (п. 3.3).

Фиг. 6.

Классификация методов оптимизации решения задачи (5).

4.1. Глобальные стохастические методы

К глобальным стохастическим методам относятся природоподобные алгоритмы, моделирующие процессы естественного отбора (генетический алгоритм [28], [35], метод дифференциальной эволюции, генетическое программирование и т.п.), а также алгоритмы, имитирующие инстинктивное поведение в природе (метод роя частиц, муравьиный алгоритм) и процессы в физике (метод имитации отжига [26]). Особенность алгоритмов состоит в решении более простой задачи, в которой используются законы природы и не учитываются особенности функционала, к которому они применяются. Такие алгоритмы естественным образом могут быть перенесены на параллельные архитектуры. Несмотря на то что они относятся к методам глобальной оптимизации, т.е. определяют область глобального экстремума, не существует теоретических оценок сходимости алгоритмов (только стохастические). Более того, настройка параметров природоподобных алгоритмов, при которых будет наблюдаться статистическая сходимость, является сложной задачей, порой сложнее самой задачи оптимизации.

4.2. Тензорное разложение

Методы покрытий нуждаются в априорной информации о функционале, экстремум которого необходимо найти. Например, основная идея метода неравномерных покрытий [32] состоит в разбиении пространства решений Q на подмножества, покрывающие пространство Q. Целевая функция на каждом подмножестве обладает свойствами (например, удовлетворение условиям липшиц-непрерывности, выпуклости, существования второй производной и т.п.), которые используются при решении задачи оптимизации, увеличивая скорость сходимости.

Тензорное разложение (ТТ) является одним из методов глобальной оптимизации и напрямую применяется к задачам максимизации параметров [37]. Он основан на свойствах TT-разложения [37], [38]

$\mathcal{T}({{i}_{1}}, \ldots ,{{i}_{d}}) = \mathop \sum \limits_{{{a}_{0}}, \ldots ,{{a}_{d}}} {{G}_{1}}({{a}_{0}},{{i}_{1}},{{a}_{1}}){{G}_{2}}({{a}_{1}},{{i}_{2}},{{a}_{2}}) \ldots {{G}_{d}}({{a}_{{d - 1}}},{{i}_{d}},{{a}_{d}}),$
и методе TT-cross аппроксимации тензора $\mathcal{T}$ [39], который строит приближение по наибольшим элементам. Здесь Gi называются ТТ-ядрами тензора $\mathcal{T}$ и являются матрицами размера ${{r}_{{i - 1}}} \times {{r}_{i}},$ ${{r}_{0}} = {{r}_{d}} = 1,$ ${{a}_{0}} = {{a}_{d}} = 1.$

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

4.2.1. Алгоритм тензорного разложения для задачи оптимизации. Задача минимизации функционала (5) может быть сведена к задаче максимизации функции $g(J(p)) = \operatorname{arcctg} J(p)$. Определим множество допустимых значений, которое могут принимать ${{p}_{k}} \in [{{a}_{k}},{{b}_{k}}]$ – компоненты вектора $p = ({{p}_{1}}, \ldots ,{{p}_{d}})$, и представим их на сетке, разбив каждый из отрезков $[{{a}_{k}},{{b}_{k}}]$ на $n$ узлов.

Пусть ${{i}_{k}}$ – произвольное значение, принимаемое параметром ${{p}_{k}}$. Взяв $g(J(p))$ от всех возможных комбинаций $p_{k}^{{{{i}_{k}}}}$, учитывая дискретность полученного ответа и представимость функционала от многих переменных в виде тензора, получаем

$\mathcal{T}({{i}_{1}}, \ldots ,{{i}_{d}}) = g(p_{1}^{{{{i}_{1}}}}, \ldots ,p_{d}^{{{{i}_{d}}}}).$

Исходная задача сводится к нахождению наибольшего элемента тензора $\mathcal{T}$. В итоге мы получим решение задачи глобальной минимизации в проекции на сетку, а итоговый результат может быть улучшен путем дальнейшей локальной минимизации.

4.3. Методы градиентного типа

Методы градиентного типа решения задачи минимизации (5) состоят в последовательном приближении решения p согласно итерационному процессу

${{p}_{{n + 1}}} = {{p}_{n}} - {{\alpha }_{n}}J{\kern 1pt} '({{p}_{n}}),\quad {{\alpha }_{n}} > 0,\quad {{p}_{0}} \in Q.$
Здесь ${{\alpha }_{n}}$ – параметр спуска, выбор которого определяет тип градиентного метода (метод наискорейшего спуска, итерации Ландвебера, минимальных ошибок, сопряженных градиентов и т.д.), $J{\kern 1pt} '({{p}_{n}}) = 2[A{\kern 1pt} '({{p}_{n}})]{\kern 1pt} {\text{*}}(A({{p}_{n}}) - f)$ – градиент целевого функционала $J(p)$. Сходимость градиентных методов к нормальному псевдорешению, а также вариации ускорения сходимости приведены в работах [20], [40]–[43] и ссылках в них. В работе [44] получено выражение для градиента $J{\kern 1pt} '$ целевого функционала:
$J{\kern 1pt} '(p) = - {{\left( {\mathop \smallint \limits_{{{t}_{0}}}^T F_{\varphi }^{{\text{T}}}(X(t),\varphi ){\Psi }(t)dt,~{{{\Psi }}_{2}}({{t}_{0}}),{{{\Psi }}_{3}}({{t}_{0}}),{{{\Psi }}_{4}}({{t}_{0}}),{{{\Psi }}_{6}}({{t}_{0}})} \right)}^{{\text{T}}}}{\kern 1pt} ,$
где ${\Psi }(t) = \left( {{{{\Psi }}_{1}}(t), \ldots ,{{{\Psi }}_{8}}(t)} \right)$ – решение сопряженной задачи
$\begin{gathered} {\dot {\Psi }} = - F_{X}^{{\text{T}}}(X(t),\varphi )\Psi (t),\quad t \in \bigcup\limits_{k = 0}^K {({{t}_{k}},{{t}_{{k + 1}}})} ,\quad {{t}_{{K + 1}}} = T; \\ \Psi (T) = 0,\quad {{[{{{\Psi }}_{i}}]}_{{t = {{t}_{k}}}}} = \frac{2}{K}({{X}_{i}}({{t}_{k}};p) - {{f}_{{ik}}}),\quad k = 1, \ldots ,K,\quad i \subset \{ 1, \ldots ,N\} . \\ \end{gathered} $
Здесь ${{F}_{\varphi }},~\;{{F}_{X}}$ – соответствующие матрицы Якоби, ${{[{{{\Psi }}_{i}}]}_{{t = {{t}_{k}}}}}$ – скачок функции ${{{\Psi }}_{i}}$ в точке tk.

5. ЗАКЛЮЧЕНИЕ

Приведены математические модели распространения нового Уханьского коронавируса COVID-2019, вспышка эпидемии которого началась в декабре 2019 г. в провинции Ухань. Модели описываются системами нелинейных обыкновенных дифференциальных уравнений, коэффициенты и начальные данные которых неизвестны или заданы усредненно. Проведен анализ идентифицируемости математической модели распространения COVID-2019 в условиях карантина, состоящей из восьми нелинейных ОДУ с 10 неизвестными коэффициентами и 4 начальными условиями, по дополнительным измерениям восприимчивых индивидуумов, инфицированных в условиях карантина и вылеченных в фиксированные моменты времени. Данная обратная задача является априорно неидентифицируемой, то есть ее решение неединственно. При фиксировании некоторых варьируемых параметров модель становится идентифицируемой. Задача идентификации параметров моделей сводится к минимизации квадратичного целевого функционала невязки. Приведены методы решения задачи минимизации, основанные на комбинации глобальных (методы покрытий, природоподобные алгоритмы, многоуровневые градиентные методы) и локальных методов (методы нулевого порядка, градиентные).

Исследуемая математическая модель схожа с математической моделью динамики туберкулеза с программами лечения [26] и коинфекции туберкулеза и ВИЧ [29]. Применение комбинированных подходов оптимизации функционала показали свою эффективность для подобного рода задач, а именно методы глобальной оптимизации определяют область глобального экстремума, а методы локальной оптимизации уточняют решения обратной задачи, “стартуя” из области глобального экстремума. Уточнение параметров в моделях развития туберкулеза в регионах Российской Федерации позволило смоделировать ситуацию на несколько лет вперед и сравнить прогноз со статистической информацией (см. фиг. 7).

Фиг. 7.

Количество чувствительных SB (бирюзовая линия), находящихся на лечении без устойчивых к лекарственным препаратам штаммов T (желтая линия) и со штаммами Tm (красная линия) в Свердловской области РФ при данных за 2010–2013 гг. с дальнейшим прогнозом на 2014–2019 гг., полученные после идентификации параметров комбинацией методов имитации отжига и градиентного спуска (тонкая линия) и тензорного разложения и градиентного спуска (толстая линия) в тыс.человек. Точками обозначены статистические данные, используемые для решения задачи идентификации параметров (красные) и для верификации алгоритмов (в цвет моделируемых величин).

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

  1. Chen T., Rui J., Wang Q., Zhao Z., Cui J., Yin L. A mathematical model for simulating the transmission of Wuhan novel Coronavirus. bioRxiv, https://doi.org/10.1101/2020.01.19.911669

  2. Huang C., Wang Y., Li X., Ren L., Zhao J. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China // Lancet (London, England). – 2020-01-24. – 24 January. – ISSN 1474-547X. https://doi.org/10.1016/S0140-6736(20)30183-5

  3. Wuhan Coronavirus (2019-nCoV) Global Cases (by JHU CSSE). Johns Hopkins University: https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6.

  4. World Health Organization. SARS (Severe Acute Respiratory Syndrome). 2019. – accessed 12th January, 2020 https://www.who.int/ith/diseases/sars/en/.

  5. Zhao S., Lin Q., Ran J., Musa S.S., Yang G., Wang W. et al. Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the earlyphase of the outbreak. Int. J. Infect. Dis. 30 Jan 2020. V. 92. P. 214–217

  6. Backer J.A., Klinkenberg D., Wallinga J. Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China, 20–28 January2020. Euro Surveill. 2020. V. 25 (5). P. 2000062.

  7. World Health Organization. Novel Coronavirus(2019-nCoV). Situation Report – 10. January 30, 2020.

  8. Ming W., Huang J., Zhang C.J.P. Breaking down of healthcare system: Mathematical modelling for controlling the novel coronavirus (2019-nCoV) outbreak in Wuhan, China. bioRxiv, https://doi.org/10.1101/2020.01.27.922443

  9. Kabanikhin S.I. Definitions and examples of inverse and ill-posed problems // J. of Inverse and Ill-Posed Problems. 2008. V. 16. P. 317–357.

  10. Waaler H.T., Geser A., Andersen S. The use of mathematical models in the study of the epidemiology of tuberculosis // Am. J. publ. Health. 1962. V. 52. P. 1002–1013.

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

  12. Perelman M.I., Marchuk G.I., Borisov S.E., Kazennukh B.Ya., Avilov K.K., Karkach A.S., Romanyukha A.A. Tuberculosis epidemiology in Russia: the mathematical model and data analysis // Russ. J. Numer. Anal. Math. Modelling. 2004. V. 19 (4). P. 305–314.

  13. Марчук Г.И. и др. Применение математического метода анализа для оценки тяжести течения, эффективности лечения, прогноза и выявления осложнений при острых пневмониях // Терапевт. арх. 1981. Т. 53. № 3. С. 3–9.

  14. Miao H., Xia X., Perelson A.S., Wu H. On identifiability of nonlinear ODE models and applications in viral dynamics // SIAM Rev. Soc. Ind. Appl. Math. 2011. V. 53 (1). P. 3–39.

  15. Kabanikhin S.I., Voronov D.A., Grodz A.A., Krivorotko O.I. Identifiability of mathematical models in medical biology // Russian Journal of Genetics: Applied Research.2016. V. 6 (8). P. 838–844.

  16. Zlojutro A., Rey D., Gardner L. Optimizing border control policies for global out-break mitigation // Scientific Reports. 2019. V. 9. P. 2216. https://rdcu.be/bniOs

  17. Tang B., Wang X., Li Q., Bragazzi N.L., Tang S., Xiao Y., Wu J. Estimation of the transmission risk of 2019-nCoV and its implication for public health interventions. SSRN: https://ssrn.com/abstract=3525558

  18. White L.F., Pagano M.A likelihood-based method for real-time estimation of the serial interval and reproductive number of an epidemic // Stat. Med. 2008. V. 27. P. 2999–3016.

  19. Chen Y., Cheng J., Jiang Y., Liu K. A time delay dynamical model for outbreak of 2019-nCoV and the parameter identification. Journal of Inverse and Ill-Posed Problems. 2020. Vol. 28, Issue 2. P. 243–250.

  20. Kaltenbacher B., Neubauer A., Scherzer O. Iterative Regularization Methods for Non-linear ill-posed Problems. Radon Series on Computational and Applied Mathematics. De Gruyter, 2008.

  21. Banks H.T., Hu Sh., Thompson W.C. Modeling and inverse problems in the presence of uncertainty. Boca Raton, London, New York, Washington, D.C.: CRC PRESS, 2014.

  22. Han X., de Vlas S.J., Fang L., Feng D., Cao W., Habbema J.D.F. Mathematical modelling of SARS and other infectious diseases in China: a review. Tropical Medicine and International Health. 2009. V. 14 (1). P. 91–100.

  23. Zhang J., Lou J., Ma Z., Wu J.A compartmental model for the analysis of SARS transmission patterns and outbreak control measures in China // Applied Mathematics and Computation. 2005. V. 162. P. 909–924.

  24. McLeod R.G., Brewster J.F., Gumel A.B., Slonowsky D.A. Sensitivity and uncertainty analyses for a SARS model with time-varying inputs and outputs. Mathematical Bio-sciences and Engineering. 2006. V. 3 (3). P. 527–544.

  25. Tahir M., Shah S.I.A., Zaman G., Khan T. Stability behavior of mathematical model MERS corona virus spread in population. Filomat. 2019. V. 33:12. P. 3947–3960.

  26. Kabanikhin S., Krivorotko O., Kashtanova V.A combined numerical algorithm for reconstructing the mathematical model for tuberculosis transmission with control programs // J. of Inverse and Ill-Posed Problems. 2018. V. 26 (1). P. 121–131.

  27. Kabanikhin S.I., Krivorotko O.I., Ermolenko D.V., Kashtanova V.N., Latyshenko V.A. Inverse problems of immunology and epidemiology // Eurasian Journal of Mathematical and Computer Applications. 2017. V. 5 (2). P. 14–35.

  28. Kabanikhin S., Krivorotko O., Takuadina A., Andornaya D., Zhang S. Geo-information system of spread of tuberculosis based on inversion and prediction. arXiv: 2002.02367 [q-bio.PE].

  29. Latyshenko V., Krivorotko O., Kabanikhin S., Zhang Sh., Kashtanova V., Yermolenko D. Identifiability analysis of inverse problems in biology // Proc. 2017 2nd International Conference on Computational Modeling, Simulation and Applied Mathematics (CMSAM). 2017. P. 567–571.

  30. Bellman R., Astrom K. On structural identifiability // Mathematical Biosciences. 1970. V. 30. № 4. P. 65–74.

  31. Bellu G., Saccomani M.P., Audoly S., D’Angiò L. DAISY: a new software tool to test global identifiability of biological and physiological systems // Comput. Methods Programs Biomed. 2007. V. 88 (1). P. 52–61.

  32. Евтушенко Ю.Г. Численный метод поиска глобального экстремума функций (перебор на неравномерной сетке) // Ж. вычисл. матем. и матем. физ. 1971. Т. 11. № 6. С. 1390–1403.

  33. Евтушенко Ю.Г., Посыпкин М.А. Метод неравномерных покрытий для решения задач многокритериальной оптимизации с гарантированной точностью // Ж. вычисл. матем. и матем. физ. 2013. Т. 53. № 2. С. 209–224.

  34. Евтушенко Ю.Г., Посыпкин М.А., Рыбак Л.А., Туркин А.В. Отыскание множеств решений систем нелинейных неравенств // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 8. С. 1248–1254.

  35. Banks H.T., Kabanikhin S.I., Krivorotko O.I., Yermolenko D.V. A numerical algorithm for constructing an individual mathematical model of HIV dynamics at cellular level // J. of Inverse and Ill-Posed Problems. 2018. V. 26 (6). P. 859–873.

  36. Zheltkova V., Zheltkov D., Grossman Z., Bocharov G.A., Tyrtyshnikov E.E. Tensor based approach to the numerical treatment of the parameter estimation problems in mathematical immunology // J. of Inverse and Ill-posed Problems. 2017. V. 26 (1). P. 51–66.

  37. Oseledets I.V., Tyrtyshnikov E.E., Breaking the curse of dimensionality, or how to use SVD in many dimensions // SIAM J. Sci. Comput. 2009. V. 31. P. 3744–3759.

  38. Oseledets I.V. Tensor-train decomposition // SIAM J. Sci. Comput. 2011. V. 33. P. 2295–2317.

  39. Tyrtyshnikov E.E. Incomplete cross approximation in the mosaic–skeleton method // Computing. 2000. V. 64. P. 367–380.

  40. Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстремальные методы решения некорректных задач. М.: Физматлит, 1988. 288 с.

  41. Васин В.В. О сходимости методов градиентного типа для нелинейных уравнений // Докл. АН. 1998. Т. 359. № 1. С. 7–9.

  42. Немировский А.С. О регуляризующих свойствах метода сопряженных градиентов на некорректных задачах// . вычисл. матем. и матем. физ. 1986. Т. 26. № 3. С. 332–347.

  43. Гасников А.В. Современные численные методы оптимизации. Метод универсального градиентного спуска. М.: МФТИ, 2018.

  44. Kabanikhin S.I., Krivorotko O.I. Identification of biological models described by systems of nonlinear differential equations // J. Inverse Ill-Posed Probl. 2015. V. 23 (5). P. 519–527.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал