\section{Модели, основанные на дифференциальных уравнениях}\label{sec_equation_models}
В данном разделе приведен обзор математических моделей, основанных на обыкновенных дифференциальных уравнениях (ОДУ), приведенных в Разделе~\ref{sec_ODE_models}, на уравнениях в частных производных (УРЧП), описанных в Разделе~\ref{sec_conv-diff}, и на стохастических дифференциальных уравнениях (СДУ), приведенных в Разделах~\ref{sec_SDE} и~\ref{sec_MFG_model}, а также описание ключевой характеристики распространения заболевания -- индекс репродукции вируса (п.~\ref{sec_reproductionNo}). В качестве примера будет представлена SEIR-HCD модель, описанная в п.~\ref{sec_R0-SEIR-HCD}.

\subsection{Введение}
Дифференциальные модели основаны на законе сохранения масс и особенностях передачи инфекции. Диаграмма~\ref{fig:hist_review} иллюстрирует развитие моделей эпидемиологии с 1760 по настоящее время.
\begin{sidewaysfigure}
    \centering
    \vspace*{12cm}
    \includegraphics[width=\textwidth]{images/Diagram_hist_review_rus.pdf}
    \caption{Диаграмма развития математических моделей эпидемиологии с 1760 года на основе камерных, имитационных моделей и их комбинации. В синих рамках приведены результаты по применению ОДУ к описанию эпидемий, в желтых рамках -- УРЧП, в оранжевых -- модели ИСП и в зеленых рамках -- АОМ.}
    \label{fig:hist_review}
\end{sidewaysfigure}

Для описания распространения новой коронавирусной инфекции, вызванной вирусом штамма SARS-CoV-2 были использованы все основные достижения, отраженные в диаграмме~\ref{fig:hist_review}. На рис.~\ref{fig:COVID_review} приведена классификация существующих математических моделей динамики COVID-19 и их взаимосвязи. Одна из простейших моделей вспышки COVID-19 описывается логистическим уравнением~(\ref{eq:logistic_eq}). Добавление состояний агентов (бессимптомные, госпитализированные, критические, умершие случаи и т.д.) к логистическому уравнению приводит к камерным SIR-моделям~(\ref{model_KermackMackendrick}) и~(\ref{eq:SEIR-HCD}). Учет пространственной неоднородности приводит к диффузионно-логистическому уравнению~(\ref{model_PDE}). Дальнейшие преобразования камерных моделей можно условно разделить на два направления: учет пространственной структуры в непрерывной (модель реакции-диффузии) и дискретной (модель конечных автоматов) постановках. Добавление управления системами состояний формирует новый класс моделей ИСП~(\ref{eq_3})-(\ref{eq_4}), а учет индивидуальных характеристик агентов -- к АОМ. Отметим, что усреднение в АОМ приводит к моделям среднего поля, а именно нелинейным цепям Маркова, в которых вероятности перехода зависят от распределения состояний агентов (подробнее см.~Раздел~\ref{sec_MFG_model}).
\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.6\textwidth]{images/Diagram_COVID_epid-models_rus.pdf}
    \caption{Взаимосвязь математических моделей распространения COVID-19 на основе камерного подхода и имитационного моделирования.}
    \label{fig:COVID_review}
\end{figure}

\subsection{SIR-модели}\label{sec_ODE_models}
Математическое моделирование в эпидемиологии началось с работы D.~Bernoulli в 1760 году, в которой была продемонстрирована эффективность вакцинации населения против ветряной оспы~\cite{Bernoulli_1760}. Впоследствии появилась серия математических моделей, основанные на законе баланса масс (см.~обзорные статьи~\cite{Bacaer_2011, Brauer_2017} и приведенную в них литературу). Работы R.~Ross в 1911~\cite{Ross_1911}, A.J.~Lotka в 1920~\cite{Lotka_1920} и V.~Volterra в 1926~\cite{Volterra_1926} (модель <<хищник-жертва>>), привели к созданию камерной SIR-модели W.O.~Kermack и A.G.~McKendrick \cite{McKendrick_1926, KermackMcKendrick_1927}
\begin{eqnarray}\label{model_KermackMackendrick}
\left\{\begin{array}{ll}
    \dfrac{dS}{dt} = -\alpha \dfrac{SI}{N}, & t> 0,\\[5pt]
    \dfrac{dI}{dt} = \alpha \dfrac{SI}{N} - \beta I,\\[5pt]
    \dfrac{dR}{dt} = \beta I,
\end{array}\right.
\end{eqnarray}
в которой популяция из $N$ особей разделена на три группы (камеры): $S$~-- восприимчивые, $I$~-- инфицированные ($I << N$) и $R$~-- вылеченные и умершие, связанные между собой вероятностными переходами $\alpha, \beta\in (0,1)$ (схема модели приведена на Рис.~\ref{fig:SEIR-HCD_scheme}а). Одним из важных результатов работы~\cite{KermackMcKendrick_1927} было введение индекса репродукции (заразности)
\[\mathcal{R}_0 = \dfrac{\alpha}{\beta},\]
которая является важнейшей характеристикой заболевания и параметром распространения эпидемии (см. подробнее раздел~\ref{sec_reproductionNo}).

\begin{sidewaysfigure}
\vspace*{12cm}
    \begin{minipage}[h]{0.26\linewidth}
    \center{\includegraphics[width=0.95\textwidth]{images/SIR_scheme.pdf}\\ а) SIR\\[5pt] \includegraphics[width=\textwidth]{images/SEIR_scheme.pdf}\\ б) SEIR} 
    \end{minipage}
    \hfill
    \begin{minipage}[h]{0.6\linewidth}
    \center{\includegraphics[width=\textwidth]{images/SEIR-HCD_scheme.pdf}\\ в) SEIR-HCD} 
    \end{minipage}
    \caption{Схемы камерных моделей (а) SIR, (б) SEIR и (в) SEIR-HCD.}
    \label{fig:SEIR-HCD_scheme}
\end{sidewaysfigure}

Отметим, что частным случаем модели~(\ref{model_KermackMackendrick}) при условии отсутствия иммунитета после перенесенного заболевания в рамках исследуемых временных промежутков (например, сезонный грипп) является SI-модель, решение которой сводится к логистическому уравнению:
\begin{eqnarray}\label{eq:logistic_eq}
    \dfrac{dI}{dt} = (\alpha-\beta)I(t) + \dfrac{\alpha}{N}I^2(t).
\end{eqnarray}

Важным свойством модели~(\ref{model_KermackMackendrick}) является выполнение закона действующих масс (закон сохранения), в рамках которого моделируемая популяция постоянна в течение всего времени $N = S(t)+I(t)+R(t)$.

В работах Е.Н.~Пелиновского и его коллег~\cite{Pelinovsky_2020, Wang_2020, Koltcova_2020} для моделирования распространения COVID-19 применяется обобщенное логистическое уравнение, описывающее рост численности заболевших. Предположение о единственности пика вспышки $I(t)$ эпидемии ограничивает применение логистической модели для описания длительного периода пандемии и учета ограничительных мер.

Для учета инкубационного периода течения COVID-19 используется модификация модели Кермака-Маккендрика SEIR-типа (схема SEIR-модели приведена на Рис.~\ref{fig:SEIR-HCD_scheme}б), которых на сегодняшний день разработано более 100 моделей~(см. например, работы~\cite{Cheng_2020, Tamm_2020, Unlu_2020, KOI_KSI_2020, Borovkov_2020, Yang_2021, Kolpakov_2021} и ссылки в них). В этих моделях популяция разделяется на группы (кроме $S$, $I$, $R$ добавляются $E$ -- бессимптомные носители, $H$ -- госпитализированные, $C$ -- критические случаи, требующие подключения аппарата искусственной вентиляции легких (ИВЛ), $D$ -- умершие в результате COVID-19, $Q$~-- помещенные на карантин и другие). Это позволяет уточнить эпидемиологическую картину в регионе за счет варьирования более детального набора коэффициентов в уравнениях. Недостатком SIR-моделей является отсутствие гибкости -- невозможность учета изменения параметров (новые мутации вируса и штамма, ограничительные меры, вакцинация). При попытке ввести в SIR-модели указанные изменения (например, сделать переменной скорость передачи инфекции $\alpha=\alpha (t)$)~\cite{Popivanov_2021}, мы сталкиваемся с неединственностью и неустойчивостью решения обратной задачи идентификации этого параметра $\alpha(t)$.

Отметим, что SIR-модели используют также для прогнозирования результатов управления развития пандемии~\cite{Silva_2021}, т.е. в правую часть уравнений добавляют кусочно-постоянную функцию управления (ограничительные меры: ношение масок, социальная дистанция, карантин). Однако и в этих случаях проблема уточнения коэффициентов моделей SIR остается открытой и требует применения методов теории обратных задач.

Также в SIR-моделях можно учитывать влияние суперраспространителей на распространение COVID-19 (инфицированные индивидуумы, у которых повышенная вирусная концентрация)~\cite{Ndairou_2020}. Однако теоретическое определение этого явления требует моделирования в масштабе отдельных людей (АОМ, см.~Раздел~\ref{sec_ABM}).

Дальнейшее развитие математических моделей можно разделить на две составляющие: введение пространственной координаты в логистические уравнения и учет дискретной пространственной неоднородности. В первом случае мы получаем новый класс математических моделей <<реакции-диффузии>> (см.~Раздел~\ref{sec_conv-diff}), а во втором -- новый подход, в котором системы дифференциальных уравнений соединяются в пространстве графовой структурой (см.~Раздел~\ref{sec_CA-ABM}).

\subsection{Модели <<реакции-диффузии>>}\label{sec_conv-diff}
В 1937 году британский ученый R.A.~Fisher~\cite{Fisher_1937} и советские математики А.Н.~Колмогоров, И.Г.~Петровский и Н.С.~Пискунов~\cite{KolmPetrPiskunov_1937} предложили и обосновали математическую модель, основанную на уравнении в частных производных параболического типа, получившую впоследствии название модели <<реакции-диффузии>>, которая была также применена при описании процессов в биологии, экологии и других приложениях:
\begin{eqnarray}\label{model_PDE}
    \dfrac{\partial u}{\partial t} = f(u) + d\dfrac{\partial^2 u}{\partial x^2}.
\end{eqnarray}
Здесь $u(x,t)$ -- вектор плотности распределения групп популяции в точке пространства $x$ и времени $t$, $d$~-- коэффициент диффузии, $f(u)$~-- функция, характеризующая характер распространения заболевания в популяции, удовлетворяющая закону сохранения масс и условиям
\[f(0)=f(1)=0,\quad f(u)>0,\;\; \text{если}\;\; 0<u<1,\]
\[f^\prime (0)>0\;\; \text{и} \;\; f^\prime (u)<f^\prime (0), \;\; \text{если}\;\; 0<u\leqslant 1.\]
Авторы работы~\cite{KolmPetrPiskunov_1937} строго доказали, что если начальное условие удовлетворяет следующим ограничениям
\[0\leqslant u(x,0)\leqslant 1,\; u(x,0)=0\;\; \forall\, x<x_1 \;\; \text{и}\;\; u(x,0)=1\;\; \forall\, x>x_2\geqslant x_1,\]
то динамика популяции по переменной $t$ описывается скоростью $v^*=2\sqrt{f^\prime(0)d}$.

Учет пространственной неоднородности позволяет более точно моделировать эпидемию от очага распространения (крупного города в стране, столицы в регионе и пр.) при известных начальных условиях. Так, в работах~\cite{Viguerie_2020, Aristov_2021, Barwolff_2021, Lau_2021_arXiv} получены оценки распространения COVID-19 в первые месяцы с начала эпидемии с учетом пассажиропотоков. Показано, что учет неоднородности влияет на характер распространения COVID-19 в крупных регионах и странах. В разделе~\ref{sec_MFG_model} будут рассмотрены стохастические дифференциальные уравнения, учитывающие пространственную неоднородность и элементы управления (задачи ИСП). Однако использование модели для описания второй и последующих волн эпидемии требует добавления уравнений в~(\ref{model_PDE}), введения множественных известных источников распространения заболевания $u(x,0)$ (обратная задача определения источника) и вычислительных ресурсов.

\subsubsection{Пространственная SIR-модель распространения инфекционного заболевания}
Приведем пример пространственной SIR-модели описания произвольного эпидемического процесса, основываясь на работах~\cite{Capasso_1997, Bratus_2010}. Введем в пространстве $\Omega$ плотности основных подгрупп популяции в момент времени $t$
\[S(t) = \int\limits_{\Omega}s(t,x)\ dx,\; I(t) = \int\limits_{\Omega}i(t,x)\ dx,\; R(t) = \int\limits_{\Omega}r(t,x)\ dx.\]

Запишем модель <<реакции-диффузии>>~(\ref{model_PDE}) SIR-типа, где описание динамики (реакция) использует соответствующие правые части SIR-моделей, а коэффициенты диффузии $d=(d_1,d_2,d_3)$ при вторых пространственных переменных полагаются постоянными величинами:
\begin{align*}
\left\{\begin{array}{ll}
    \dfrac{\partial s}{\partial t} = -\alpha \dfrac{si}{N} + d_1 s_{xx}, & t> 0,\, x\in\Omega,\\[5pt]
    \dfrac{\partial i}{\partial t} = \alpha \dfrac{si}{N} - \beta i + d_2 i_{xx},\\[5pt]
    \dfrac{\partial r}{\partial t} = \beta i + d_3 r_{xx},
\end{array}\right.
\end{align*}

\subsection{Стохастические модели}\label{sec_SDE}
Стохастические модели эпидемии изучались в различных сетях~\cite{May_Lloyd_2001, Pastor_2001, Newman_2002}, а также применялись к моделированию пандемии COVID-19~\cite{Wong_Collins_2020}.

В данном разделе пойдет речь о трех различных методах построения стохастических моделей эпидемий, которые непосредственно связаны с их детерминированными аналогами~\cite{Allen2008}: 1) дискретная временная марковская цепная модель (Раздел~\ref{sec_DTMC}), 2) цепная модель Маркова с непрерывным временем (Раздел~\ref{sec_CTMC}), и 3) модель, основанная на системе СДУ (Раздел~\ref{sec_SDE_system}). Эти стохастические процессы различаются по исходным предположениям относительно времени и состояния переменных. В модели 1 время и переменные состояния дискретны. В модели 2 время является непрерывным, а переменная состояния -- дискретной. Наконец, модель СДУ основана на диффузионном процессе, где и время, и переменные состояния являются непрерывными. Описание и свойства стохастических моделей эпидемий приведем для SIR-модели~(\ref{model_KermackMackendrick}).

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

Всюду далее каллиграфические символы $\mathcal{S}(t)$, $\mathcal{I}(t)$, $\mathcal{R}(t)$ будут означать случайные величины (за исключением базового индекса репродукции $\mathcal{R}_0$).

\subsubsection{Дискретная временная марковская цепная модель}\label{sec_DTMC}
Пусть $\mathcal{S}(t)$, $\mathcal{I}(t)$ и $\mathcal{R}(t)$~-- дискретные случайные величины, описывающие число восприимчивых, инфицированных и умерших/вылеченных индивидуумов в момент времени $t\in \{0, \Delta t, 2\Delta t, \ldots\}$, соответственно, т.е. $\mathcal{S}(t), \mathcal{I}(t), \mathcal{R}(t) \in \{0,1,2,\ldots,N\}$.

В данном примере мы имеем дело с двумя независимыми случайными величинами $\mathcal{S}(t)$ и $\mathcal{I}(t)$. Случайная величина $\mathcal{R}(t) = N - \mathcal{S}(t) - \mathcal{I}(t)$. Случайный двумерный процесс $\{\mathcal{S}(t), \mathcal{I}(t)\}_{t=0}^\infty$ характеризуется совместной функцией вероятности
\[p_{(s,i)}(t) = \mbox{Prob}\ \{\mathcal{S}(t) = s, \mathcal{I}(t) = i\}.\]
Данный нелинейный процесс обладает свойством <<марковости>> и однородности по времени.

\begin{definition}
\textit{Свойство <<марковости>>} стохастического процесса состоит в том, что процесс в момент времени $t+\Delta t$ зависит только он своего состояния в предыдущий момент времени $t$.
\end{definition}

\begin{definition}
Если вероятность перехода $p_{ji}(t+\Delta t, t) = \mbox{Prob}\ \{\mathcal{I}(t+\Delta t) = j\ | \ \mathcal{I}(t)=i\}$ из состояния $i$ в состояние $j$ не зависит от $t$, т.е. $p_{ji}(\Delta t)$, то случайный процесс $\{\mathcal{I}(t)\}_{t=0}^\infty$ является \textit{однородным по времени}.
\end{definition}

Вероятности перехода между состояниями основываются на предположениях SIR-модели~(\ref{model_KermackMackendrick}). Предположим, что $\Delta t$ достаточно мал, чтобы за это время произошло более одного события (заражение, выздоровление, смерть). Тогда вероятность перехода между состояниями имеет вид:
\[p_{(s+k, i+j), (s,i)}(\Delta t) = \mbox{Prob}\{(\Delta \mathcal{S}, \Delta\mathcal{I}) = (k,j) \ | \ (\mathcal{S}(t), \mathcal{I}(t)) = (s,i)\}.\]
Здесь $\Delta\mathcal{S} = \mathcal{S}(t+\Delta t) - \mathcal{S}(t)$. Пользуясь формулировкой SIR-модели~(\ref{model_KermackMackendrick}), получим
\begin{align*}
p_{(s+k, i+j), (s,i)}(\Delta t) = \left\{\begin{array}{ll}
\dfrac{\alpha is}{N}\Delta t, & (k,j) = (-1,1),\\
\beta i\Delta t, & (k,j) = (0,-1),\\
1 - \dfrac{\alpha is}{N}\Delta t - \beta i\Delta t, & (k,j) = (0,0),\\
0, & \text{иначе}.
\end{array}\right.
\end{align*}
Таким образом, вероятность инфицирования (переход $(s,i)\to (s-1, i+1)$) равна $\dfrac{\alpha is}{N}\Delta t$, а вероятность выздоровления или смерти (переход $(s,i)\to (s, i-1)$) равна $\beta i\Delta t$. Наконец, вероятность отсутствия изменений (переход $(s,i)\to (s, i)$) имеет вид $1 - \dfrac{\alpha is}{N}\Delta t - \beta i\Delta t$. Т.к. шаг по времени $\Delta t$ достаточно мал, каждая вероятность перехода между состояниями лежит в интервале $[0,1]$.

Применяя свойство <<марковости>> случайного процесса $\{\mathcal{S}(t), \mathcal{I}(t)\}_{t=0}^\infty$, получим дифференциальное уравнение для вероятности $p_{(s,i)}(t+\Delta t)$:
\begin{align}\label{eq_DTMC_model}
\begin{array}{ll}
	p_{(s,i)}(t+\Delta t) & = p_{(s+1,i-1)}(t) \dfrac{\alpha}{N}(i-1)(s+1)\Delta t + p_{(s,i+1)}(t) \beta (i+1)\Delta t + \\[5pt]
	& + p_{(s,i)} (t) \left( 1 - \left[\dfrac{\alpha}{N}is + \beta i\right]\Delta t\right).
\end{array}
\end{align}

Несложно показать, что состояние без заболевания $(N,0)$ является поглощающим, т.е. все траектории процесса асимптотически стремятся к $(N,0)$.

\subsubsection{Цепная модель Маркова с непрерывным временем}\label{sec_CTMC}
Пусть переменные состояния $\mathcal{S}(t)$, $\mathcal{I}(t)$ и $\mathcal{R}(t)$ по-прежнему являются дискретными случайными величинами, т.е. $\mathcal{S}(t), \mathcal{I}(t), \mathcal{R}(t) \in \{0,1,2,\ldots,N\}$, а время $t\in (0,T)$ непрерывно.

Вероятность перехода между состояниями процесса $\{\mathcal{S}(t), \mathcal{I}(t)\}_{t\in (0,T)}$ характеризуется следующей системой, как и в предыдущем случае, с добавлением $o(\Delta t)$:
\begin{align*}
p_{(s+k, i+j), (s,i)}(\Delta t) = \left\{\begin{array}{ll}
\dfrac{\alpha is}{N}\Delta t + o(\Delta t), & (k,j) = (-1,1),\\
\beta i\Delta t + o(\Delta t), & (k,j) = (0,-1),\\
1 - \dfrac{\alpha is}{N}\Delta t - \beta i\Delta t + o(\Delta t), & (k,j) = (0,0),\\
0, & \text{иначе}.
\end{array}\right.
\end{align*}
В результате мы получим схожее уравнение с~(\ref{eq_DTMC_model}) с добавлением $o(\Delta t)$ в правую часть. Полагая $\Delta t \to 0$ в итоговом уравнении для переходной вероятности, получим прямое уравнение Колмогорова:
\begin{align*}
\begin{array}{ll}
	\dfrac{dp_{(s,i)}}{dt} & = p_{(s+1,i-1)} \dfrac{\alpha}{N}(i-1)(s+1) + p_{(s,i+1)} \beta (i+1) - p_{(s,i)} \left( \dfrac{\alpha}{N}is + \beta i\right).
\end{array}
\end{align*}

В работе 2011 года~\cite{Sazonov_2011} показано, что время достижения вспышки эпидемии не сильно зависит от начальной траектории, описываемой в рамках SI марковской модели. В работе построена гибридная модель, состоящая из стохастического марковского процесса на начальном этапе моделирования с целью получения распределения исходных условий распределения S и I популяций, и детерминированной SIR модели для описания последующей эволюции системы. В работе тех же авторов 2016 года~\cite{Sazonov_2016} построена модель динамики эпидемии в стохастических взаимодействующих центрах популяции в сочетании с случайной миграцией на основе марковских цепей SIR типа. Получено среднее приближение поля для общего процесса миграции и предложен приближенный метод, позволяющий вычислять статистические моменты для сетей с густонаселенными центрами. Рассматриваемые модели решаются численно на основе метода Монте Карло, что увеличивает их вычислительную сложность. В работе 2021 года~\cite{Sazonov_2021_HIV_num} предложен гибридный метод численного решения задачи идентификации параметров кинетики жизненного цикла ВИЧ на клеточном уровне, в основе которого лежит марковский процесс и анализ чувствительности стохастической модели.

\subsubsection{Стохастические дифференциальные уравнения}\label{sec_SDE_system}
Пусть состояния SIR-системы и время непрерывны, т.е. $\mathcal{S}(t), \mathcal{I}(t), \mathcal{R}(t)\in [0,N]$ и $t\in (0,T)$.
Предположим, что независимые случайные величины $\{(\mathcal{S}(t),\mathcal{I}(t))\}$ имеют нормальное распределение (или асимптотически стремятся к такому).

Обозначим $\Delta X(t) = \left( \Delta\mathcal{S}, \Delta\mathcal{I}\right)^T$. Тогда математическое ожидание $\Delta X(t)$ порядка $\Delta t$ записывается в виде
\begin{align*}
	\mathbb{E}(\Delta X(t)) = \left(\begin{array}{c}
	-\dfrac{\alpha}{N}\mathcal{S}\mathcal{I} \\
	\dfrac{\alpha}{N}\mathcal{S}\mathcal{I} - \beta\mathcal{I}
	\end{array}\right) \Delta t.
\end{align*}

Матрица ковариации для $\Delta X(t)$ может быть записана в первом приближении:
\[V(\Delta X(t)) = \mathbb{E}(\Delta X(t)[\Delta X(t)]^T) - \mathbb{E}(\Delta X(t))\mathbb{E}(\Delta X(t))^T \approx \mathbb{E}(\Delta X(t)[\Delta X(t)]^T)\]
или
\begin{align*}
	V(\Delta X(t)) = \left(\begin{array}{cc}
	\dfrac{\alpha}{N}\mathcal{S}\mathcal{I} & -\dfrac{\alpha}{N}\mathcal{S}\mathcal{I} \\
	-\dfrac{\alpha}{N}\mathcal{S}\mathcal{I} & \dfrac{\alpha}{N}\mathcal{S}\mathcal{I} + \beta\mathcal{I}
	\end{array}\right) \Delta t.
\end{align*}

В силу того, что случайная величина $\Delta X(t)$ имеет нормальное распределение, то имеет место следующее разложение
\begin{align}\label{eq_vector_SDE}
	X(t+\Delta t) = X(t) + \Delta X(t) \approx X(t) + \mathbb{E}(\Delta X(t)) + \sqrt{V(\Delta X(t))}.
\end{align}
В силу того, что $V>0$ и симметричная, существует однозначный корень $\sqrt{V} = B\sqrt{\Delta t}$. Система уравнений~(\ref{eq_vector_SDE}) является аппроксимацией Эйлера системы СДУ в смысле Ито. При достаточно гладких коэффициентах решение $X(t)$ системы~(\ref{eq_vector_SDE}) сходится к решению системы СДУ:
\begin{align}\label{eq_SDE_system}
	\left\{\begin{array}{ll}
		d\mathcal{S}= -\dfrac{\alpha}{N}\mathcal{S}\mathcal{I}dt + B_{11}dW_1 + B_{12}dW_2,\\[5pt]
		d\mathcal{I} = \left(\dfrac{\alpha}{N}\mathcal{S}\mathcal{I} - \beta\mathcal{I}\right) dt + B_{21}dW_1 + B_{22}dW_2,
	\end{array}\right.
\end{align}
где $W_1$ и $W_2$ -- два независимых винеровских процесса, $B=(B_{ij})$.

\begin{definition}
	Случайный процесс $W_t$, $t\geq 0$, называется винеровским процессом, если
	\begin{enumerate}
	\itemsep0em
	    \item $W_0 = 0$ почти достоверно.
	    \item $W_t$ -- процесс с независимыми приращениями.
	    \item $W_{t+\Delta t} - W_t \sim \mathbf{N}(0,\Delta t)$. Здесь $\mathbf{N}(0,\Delta t)$~-- нормальное распределение со средним 0 и дисперсией $\Delta t$.
	\end{enumerate}
\end{definition}

На рисунке~\ref{fig:sir_determined_and_stochastic} приведены решения уравнений для SIR модели в детерминированной постановке~(\ref{model_KermackMackendrick}) (отмечены черной линией) и в стохастической~(\ref{eq_SDE_system}) для восприимчивых $S(t)$ (рис.~\ref{fig:sir_determined_and_stochastic}а) и инфицированных $I(t)$ (рис.~\ref{fig:sir_determined_and_stochastic}б) индивидуумов в случае размера моделируемой популяции $N=100$ человек, $\alpha = 1$, $\beta = 0.5$, $S_0 = 98$, $I_0 = 2$, шаг по времени брался равномерным $\Delta t = 0.01$. Для решения системы~(\ref{model_KermackMackendrick}) применялась схема Эйлера 1-го порядка аппроксимации, решение системы~(\ref{eq_SDE_system}) строилось на основе схемы Эйлера-Маруямы:
\begin{align*}
	S_{i+1} & = S_i - \dfrac{\alpha}{N}S_i I_i\Delta t + B_{11}\sqrt{\Delta t}\,\mbox{rand}_1 + B_{12}\sqrt{\Delta t}\,\mbox{rand}_2,,\\[5pt]
	I_{i+1} & = I_i + (\dfrac{\alpha}{N}S_i I_i - \beta I_i)\Delta t + B_{21}\sqrt{\Delta t}\,\mbox{rand}_1 + B_{22}\sqrt{\Delta t}\,\mbox{rand}_2,
\end{align*}
где $\mbox{rand}_1\sim \mathbf{N}(0,1)$ и $\mbox{rand}_2\sim \mathbf{N}(0,1)$~-- две независимые случайные величины.

\begin{figure}[!ht]
\centering
\begin{minipage}{0.9\linewidth}
\centering{\includegraphics[width=\linewidth]{images/SIR_S.png} а) Количество восприимчивых индивидуумов $S(t)$}
\end{minipage}
\vfill
\begin{minipage}{0.9\linewidth}
\centering{\includegraphics[width=\linewidth]{images/SIR_I.png} б) Количество инфицированных индивидуумов $I(t)$}
\end{minipage}
\caption{Моделирование изменения $S(t)$ и $I(t)$ в детерминированной~(\ref{model_KermackMackendrick}) (черная линия) и стохастической~(\ref{eq_SDE_system}) (цветные линии) SIR модели при условиях: $N=100$, $\alpha = 1$, $\beta = 0.5$, $S_0 = 98$, $I_0 = 2$, $\Delta t = 0.01$.}
\label{fig:sir_determined_and_stochastic}
\end{figure}

Траектории решения системы~(\ref{eq_SDE_system}) достаточно близки к решению детерминированной SIR модели~(\ref{model_KermackMackendrick}).

В работе~\cite{Allen2008} показано, что в случае достаточно большой популяции $N\geqslant 100$, индекса репродукции вируса $\mathcal{R}_0$ и если вспышка начинается с достаточным количеством инфицированных особей, то результаты для стохастической SIR эпидемии находятся в близком согласии с предсказаниями детерминированной эпидемической SIR модели; болезнь становится эндемичной.

\subsubsection{Случайные процессы в эпидемиологии}
Руководствуясь случайностью изучаемых процессов, анализ эпидемии возможен и без явного построения стохастических моделей. Известно, что плотность потока частиц в размножающей среде при достаточно широких условиях асимптотически экспоненциальна по времени $t$ с некоторым параметром $\lambda$, т.е. с показателем $\lambda t$. В работах~\cite{Lotova_2020, Lotova_2021} показано, что если среда случайна, то параметр $\lambda$ -- случайная величина, и для оценки временной асимптотики среднего (по реализациям среды) числа частиц можно в некотором приближении усреднять экспоненту по распределению $\lambda$. В предположении гауссовости этого распределения таким образом получается асимптотическая “сверхэкспоненциальная” оценка среднего потока, выражаемая экспонентой с показателем $t\mathbb{E}\lambda + t^2D\lambda/2$. Для численной экспериментальной проверки такой оценки разработано вычисление вероятностных моментов случайного параметра $\lambda$ на основе рандомизации фурье-приближений специальных нелинейных функционалов.

Отмечено, что результаты в статье имеют широкое приложение. В частности, согласно статистике Всемирной организации здравоохранения (ВОЗ), сверхэкспоненциальное поведение проявляла пандемия COVID-19, развивающаяся во всем мире. А именно, количество наблюдений (по дням) аппроксимируется с точностью до 2\% с 9 марта 2020 года по 21 марта 2020 года.

\subsection{Индекс репродукции вируса $\mathcal{R}_0$ для дифференциальных моделей}\label{sec_reproductionNo}
Общее определение индекса репродукции (basic reproduction number)~\cite{R0_wiki}: индекс репродукции вируса $\mathcal{R}_0$ определяется как среднее количество людей, которых заражает активный инфицированный, попавший в полностью неиммунизированное окружение при отсутствии специальных эпидемиологических мер, направленных на предотвращение распространения заболевания.

В разделе~\ref{sec_R0} мы покажем, что индекс репродукции $\mathcal{R}_0$ является границей устойчивости состояния равновесия SIR-системы при отсутствии инфицированных. Если $\mathcal{R}_0>1$, то на начальном этапе число заболевших будет расти экспоненциально. Если $\mathcal{R}_0 \in (0,1)$, то небольшое количество инфицированных людей, попавших в полностью восприимчивую популяцию, в среднем не смогут сохранить свою группу, и эпидемии не будет.

\subsubsection{Алгоритм вычисления индекса репродукции}\label{sec_R0}
%We focuses on compartmental models for disease transmission based on ordinary differential equations (ODE) and next-generation method proposed by van den Driessche and Watmough~\cite{Driessche_2008}. A compartment is called a \textit{disease} compartment if the individuals therein are infected. Note that this use of the term disease is broader than the clinical definition and includes asymptomatic stages of infection as well as symptomatic. Based on approach from~\cite{Driessche_2008} we consider $n$ disease compartments and $m$ nondisease compartments, and let $x\in \mathbb{R}^n$ and $y\in\mathbb{R}^m$ be the subpopulations in each of these compartments. Then the mathematical model of disease transmission based on ODE has the form
Опишем метод вычисления индекса репродукции, предложенный van den Driessche и Watmough~\cite{Driessche_2008} для SIR-моделей. Разделим всю популяцию на две категории: $n$ групп инфицированных и $m$ групп не инфицированных. Векторы $x\in \mathbb{R}^n$ и $y\in\mathbb{R}^m$ определяют количество индивидуумов в каждой из двух категорий, например, в модели SEIR-HCD $x=(E,I,H,C)\in \mathbb{R}^4$, $y=(S,R,D)\in \mathbb{R}^3$ (см. Раздел~\ref{sec_R0-SEIR-HCD}). Тогда SIR-модель распространения инфекционного заболевания принимает вид
\begin{eqnarray}\label{eq:ODE_2compartment}
\begin{array}{ll}
    \dot{x}_i = \mathcal{F}_i(x,y) - \mathcal{V}_i(x,y),\quad &  i=1,\ldots,n,\\
    \dot{y}_j = g_j(x,y), & j=1,\ldots,m.
\end{array}    
\end{eqnarray}
Здесь $\mathcal{F}_i$ -- скорость изменения количества инфицированных и $\mathcal{V}_i$ -- скорость изменения умерших, выздоровевших и тех, кто заболевает в $i$-й группе.

%The derivation of the basic reproduction number is based on the linearization of the ODE model about a disease-free equilibrium $(0,y_0)$. In paper~\cite{Driessche_2008} five assumptions are proposed to ensure the existence of this equilibrium and to ensure the model~(\ref{eq:ODE_2compartment}) is well posed. Then we construct matrices
Вывод индекса репродукции заключается в линеаризации системы~(\ref{eq:ODE_2compartment}) в окрестности состояния равновесия в случае отсутствия инфекции $(0,y_0)$. В работе~\cite{Driessche_2008} были введены пять условий существования этого состояния равновесия, согласованных с законом сохранения масс для модели модели~(\ref{eq:ODE_2compartment}). Вычислим собственные значения матрицы $FV^{-1}$, где
\[F_{ij} = \dfrac{\partial \mathcal{F}_i}{\partial x_j}(0,y_0),\quad V_{ij} = \dfrac{\partial \mathcal{V}_i}{\partial x_j}(0,y_0).\]
%And the reproduction number $\mathcal{R}_0 = \max \lambda_i(FV^{-1})$ is a nonnegative~\cite{Berman_1970}, where $\lambda(FV^{-1})$ is the eigenvalue of matrix $FV^{-1}$. There is a nonnegative eigenvector $\omega$ associated with $\mathcal{R}_0$. This eigenvector is in some sense the distribution of infected individuals that produces the greatest number, $\mathcal{R}_0$, of secondary infections per generation.
Тогда индекс репродукции $\mathcal{R}_0 = \max \lambda_i(FV^{-1})$ является неотрицательным и соответствующий ему собственный вектор $\omega$ состоит из неотрицательных компонент~\cite{Berman_1970}. Компоненты вектора $\omega$ можно интерпретировать как распределение инфицированных индивидуумов, вызывающих наибольшее количество $\mathcal{R}_0$ вторичных инфекций в поколении.

%In~\cite{Driessche_2008} the local stability of disease-free equilibrium theorem is proved, i.e. the disease-free equilibrium $(0,y_0)$ is locally asymptotically stable if $\mathcal{R}_0 < 1$, but unstable if $\mathcal{R}_0 > 1$. 
В~\cite{Driessche_2008} доказана теорема о локальной устойчивости неинфицированного состояния равновесия системы~(\ref{eq:ODE_2compartment}), а именно состояние равновесия $(0, y_0)$ локально асимптотически устойчиво, если $\mathcal{R}_0<1$, и неустойчиво, если $\mathcal{R}_0>1$. В следующем разделе приведен вывод индекса репродукции для SEIR-HCD модели распространения COVID-19.

\subsubsection{Индекс репродукции вируса для модели SEIR-HCD}\label{sec_R0-SEIR-HCD}
Используя описанный алгоритм, выведем индекс репродукции для модели SEIR-HCD распространения COVID-19~\cite{Unlu_2020, KOI_KSI_2020}
\begin{eqnarray}\label{eq:SEIR-HCD}
    \left\{\begin{array}{ll}
        \dfrac{dS}{dt} = -\tilde{a}(t)\left( \dfrac{\alpha_I(t) S(t)I(t)}{N} + \dfrac{\alpha_E(t) S(t)E(t)}{N}\right) + \dfrac{1}{t_{imm}}R(t), \\
        \dfrac{dE}{dt} = \tilde{a}(t)\left( \dfrac{\alpha_I(t) S(t)I(t)}{N} + \dfrac{\alpha_E(t) S(t)E(t)}{N}\right) - \dfrac{1}{t_{inc}}E(t), \\
        \dfrac{dI}{dt} = \dfrac{1}{t_{inc}}E(t) - \dfrac{1}{t_{inf}}I(t), \\
        \dfrac{dR}{dt} = \dfrac{\beta (t)}{t_{inf}}I(t) + \dfrac{1-\varepsilon_{HC}(t)}{t_{hosp}}H(t) - \dfrac{1}{t_{imm}}R(t), \\
        \dfrac{dH}{dt} = \dfrac{1-\beta (t)}{t_{inf}}I(t) + \dfrac{1-\mu (t)}{t_{crit}}C(t) - \dfrac{1}{t_{hosp}}H(t), \\
        \dfrac{dC}{dt} =  \dfrac{\varepsilon_{HC}(t)}{t_{hosp}}H(t) - \dfrac{1}{t_{crit}}C(t), \\
        \dfrac{dD}{dt} =  \dfrac{\mu (t)}{t_{crit}}C(t)
    \end{array}\right.
\end{eqnarray}
с начальными условиями
\begin{equation}\label{eq:SEIR-HCD_init_cond}
\begin{split}
    S(t_0)=N-E_0-I_0-R_0-H_0-C_0-D_0,\; E(t_0)=E_0,\\
    I(t_0)=I_0,\; R(t_0)=R_0,\; H(t_0)=H_0,\; C(t_0)=C_0,\; D(t_0)=D_0.
\end{split}
\end{equation}

Схема модели~(\ref{eq:SEIR-HCD}) приведена на Рис.~\ref{fig:SEIR-HCD_scheme}в, а описание и значения параметров и начальных условий для Новосибирской области приведены в таблице~\ref{tab_parameters} (начальный момент времени полагается 15.04.2020). В SEIR-HCD модели происходит перемещение бессимптомной популяции $E(t)$ после $t_{inc}$ дней в симптоматическую $I(t)$. Инфицированные индивидуумы после $t_{inf}$ дней выздоравливают с вероятностью $\beta$ и госпитализируются $H(t)$ с вероятностью $1-\beta$. Затем госпитализированные могут выздоравливать или нуждаться в подключении аппарата ИВЛ $C(t)$. В модели только критические случаи могут умереть $D(t)$ с вероятностью $\mu$.

Разделяя всю популяцию $N$ в модели~(\ref{eq:SEIR-HCD}) на инфицированных $x=(E,I,H,C)\in \mathbb{R}^4$ и неинфицированных $y=(S,R,D)\in \mathbb{R}^3$, получим следующие вектор-функции, согласно~(\ref{eq:ODE_2compartment}):
\begin{eqnarray*}
\begin{array}{c}
\mathcal{F}=\left(\begin{array}{c}
    \tilde{a}(t) \left(\frac{\alpha_I(t) S(t)I(t)}{N} + \frac{\alpha_E(t) S(t)E(t)}{N}\right)\\
    0 \\ 0 \\ 0
\end{array}\right), \\
\mathcal{V}=\left(\begin{array}{c}
    \frac{1}{t_{inc}}E(t)\\
    -\frac{1}{t_{inc}}E(t) + \frac{1}{t_{inf}}I(t) \\
    -\frac{1-\beta}{t_{inf}}I(t) - \frac{1-\mu}{t_{crit}}C(t) + \frac{1}{t_{hosp}}H(t) \\
    -\frac{\varepsilon_{HC}}{t_{hosp}}H(t) + \frac{1}{t_{crit}}C(t)
\end{array}\right),
\end{array}
\end{eqnarray*}
где $\tilde{a}(t) = \dfrac{5-a(t-\tau)}{5}$.

\begin{center}
\begin{table}[h!]
\caption{Описание и значения параметров SEIR-HCD модели}\label{tab_parameters}
 \begin{tabular}{|p{1.7cm}|p{8cm}|p{1.7cm}| } 
 \hline
 Параметр & Описание & Значение\\ 
  \hline
$a(t)$ & Индекс самоизоляции от Яндекса & $(0, 5)$ \\ 
$\alpha_E(t)$ & Параметр заражения между бессимптомной и восприимчивой группами населения ($\alpha_E >> \alpha_I$) & $(0, 1)$ \\
$\alpha_I(t)$ & Параметр заражения между инфицированным и восприимчивым населением & $(0,1)$\\
$\beta$ & Доля инфицированных, которая переносит заболевание без осложнений & $(0,1)$\\
$\varepsilon_{HC}$ & Доля госпитализированных случаев, которым требуется подключение ИВЛ & $(0,1)$\\
$\mu$ & Доля смертельных случаев в результате COVID-19 & $(0,0.5)$\\
$\tau$ & Латентный период (характеризует запаздывание наступления заразности) & 2 дня \\
$t_{inc}$ & Длительность инкубационного периода & 2-14 дня \\
$t_{inf}$ & Длительность периода инфицирования & 2.5-14 дня \\
$t_{hosp}$ & Длительность периода госпитализации & 4-5 дня \\
$t_{crit}$ & Длительность использования аппарата ИВЛ & 10-20 дня \\
$t_{imm}$ & Длительность гуморального иммунитета & 120-180 дней \\
$N$ & Население в Новосибирской области (человек) & 2798170 \\
$E_0$ & Начальное количество бессимптомных носителей & (1, 5000) \\
$I_0$ & Начальное количество инфицированных случаев & (1, 5000) \\
$R_0$ & Начальное количество вылеченных случаев & (1, 100) \\
$H_0$ & Начальное количество госпитализированных & 133 \\
$C_0$ & Начальное количество критических случаев & 1 \\
$D_0$ & Начальное количество смертей & 1 \\
\hline
\end{tabular}
\label{interventions}
\end{table}
\end{center}

Состояние равновесия системы~(\ref{eq:SEIR-HCD})-(\ref{eq:SEIR-HCD_init_cond}) в случае отсутствия инфицированных индивидуумов есть $\mathcal{E}_0=(N,0,0,0,0,0,0)$. Тогда матрицы $F$ и $V$ имеют вид:
\begin{eqnarray*}
F=\left(\begin{array}{cccc}
    \tilde{a}\alpha_E &  \tilde{a}\alpha_I & 0 & 0\\
    0 & 0 & 0 & 0\\
    0 & 0 & 0 & 0\\
    0 & 0 & 0 & 0
\end{array}\right), \quad
V=\left(\begin{array}{cccc}
    \frac{1}{t_{inc}} & 0 & 0 & 0\\
    -\frac{1}{t_{inc}} & \frac{1}{t_{inf}} & 0 & 0\\
    0 & -\frac{1-\beta}{t_{inf}} & \frac{1}{t_{hosp}} & -\frac{1-\mu}{t_{crit}}\\
    0 & 0 & -\frac{\varepsilon_{HC}}{t_{hosp}} & \frac{1}{t_{crit}}
\end{array}\right).
\end{eqnarray*}

Вычисляя максимальное собственное значение матрицы $FV^{-1}$, получим индекс репродукции для SEIR-HCD модели
\begin{eqnarray}\label{eq:R0_SEIR-HCD}
    \mathcal{R}_0(t) = \tilde{a}(t) \left( \alpha_E(t) t_{inc} + \alpha_I(t) t_{inf}\right),
\end{eqnarray}
первое слагаемое в котором характеризует заразность от бессимптомной части популяции, а второе слагаемое -- от симптомной.

\begin{lemma}\label{lemma_R0_Lyapunov_stability}
     В случае постоянных параметров модели (\ref{eq:SEIR-HCD}) и начальных условий~(\ref{eq:SEIR-HCD_init_cond}) состояние равновесия $\mathcal{E}_0$ системы неустойчиво.
     
     При отсутствии слагаемого $\dfrac{R(t)}{t_{imm}}$, характеризующего потерю иммунитета после выздоровления, состояние равновесия $\mathcal{E}_0$ устойчиво по Ляпунову при $\mathcal{R}_0<1$, определяемым формулой~(\ref{eq:R0_SEIR-HCD}).
\end{lemma}

\subsection{Выводы}
Преимущество использования дифференциальных моделей для описания распространения эпидемий (в том числе COVID-19) состоит в учете особенностей инфекции (наличие инкубационного периода, взаимосвязь инфицированных и критических случаев и т.д.) в параметрах моделей и закона сохранения масс (размера популяции), а в случае добавления переменной $x$ (модели <<реакции-диффузии>>) -- учет пространственной неоднородности. Такие модели качественно описывают вспышку эпидемии, но не обладают достаточной гибкостью. Учет изменения параметров для описания ограничительных мер, новых мутаций вируса приводит к неединственности и неустойчивости решения задачи их идентификации.

В стохастических SIR моделях эндемическое равновесие не достигается, поскольку существует конечная вероятность того, что число инфицированных особей в системе упадет ниже единицы. Таким образом, в данной системе патоген может не размножаться. В случае достаточно большой популяции, индекса репродукции вируса $\mathcal{R}_0>1$ и если вспышка начинается с достаточным количеством инфицированных особей, то результаты для стохастической SIR эпидемии находятся в близком согласии с предсказаниями детерминированной эпидемической SIR модели; болезнь становится эндемичной.